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Final  Technical  Report 


A.  Objectives  of  the  research  effort 

As  discussed  in  earlier  Reports,  our  research  effort  was  motivated  by  a  set  of 
experiments  by  Mark  Debe  and  coworkers,  which  explored  various  features  of  the  growth 
of  crystalline  thin  films  of  Copper  Phthalocyanine  (CuPc)  via  physical  vapor  transport 
(PVT)  using  a  closed-cell  ampoule  [1. 1-1.4].  Films  grown  under  a  variety  of  conditions  in 
microgravity  on  the  Space  Shuttle  in  low  Earth  orbit  (LEO)  were  compared  to  ground 
control  experiments  at  unit  gravity.  The  major  finding  is  that  the  center  region  of  two  films 
of  CuPc  produced  in  LEO  are  strikingly  different  from  the  edges  and  from  both  the  centers 
and  the  edge  region  of  the  13  ground-grown  control  films  of  equivalent  thickness.  In 
particular,  the  crystalline  density  is  at  least  a  factor  of  two  greater  at  the  center  than  at  the 
edges  or  for  the  ground  grown  controls.  The  LEO  central  regions  consisted  of  well- 
oriented,  densely  packed  columnar  structures-a  morphology  with  desirable  optical 
properties— while  the  unit  gravity  samples  included  only  widely  spaced  narrow  whiskers 
that  appeared  to  be  tangled  and  not  weU-oriented.  In  addition  there  is  some  evidence  that  the 
imderlying  crystalline  polymorph  may  also  vary  between  the  two  sets  of  samples. 

TTie  objective  of  the  present  research  project  was  to  imderstand  the  nature  of  the 
morphological  selection  in  the  PVT  crystal  growth  process.  One  major  goal  has  been  to 
assess  from  a  theoretical  perspective  what  role  gravity  could  play  in  producing  the  different 
morphologies  seen  in  LEO.  A  second  major  goal  has  been  to  determine  how  the  very 
anisotropic  (platelike)  shape  of  the  CuPc  molecule  influences  the  crystal  growth 
morphologies  both  on  earth  and  in  space.  As  indicated  in  detail  below  we  have  made  good 
progress  towards  both  these  goals,  though  more  work  needs  to  be  done.  Because  of  the 
third  year  of  our  contract  was  not  funded,  we  were  not  able  to  pursue  some  promising 
leads,  and  must  leave  some  work  to  the  future. 

The  most  obvious  way  gravity  could  play  a  role  is  through  its  effect  on  convection  in 
the  cell.  Accordingly,  we  devoted  considerable  effort  to  determining  whether  the  geometry 
of  the  PVT  cell,  thermal  gradients,  and  the  density  of  enclosed  gases  are  sufficient  to 
induce  convection  at  unit  gravity.  We  have  concluded  that  it  is  unlikely  that  significant 
convection  can  arise  as  a  result  of  the  one-dimensional  (vertical)  temperature  gradient  in 
Debe's  experimental  cell.  A  summary  of  our  conclusions  is  given  in  Section  B.l,  and  a 
detailed  discussion  is  attached  as  Appendix  A  of  this  Report 

A  major  problem  we  faced  in  t^ing  to  assess  the  role  of  gravity  was  that  very  little  is 
known  about  many  aspects  of  thin  film  growth  of  large  organic  molecules  even  on  earth. 
Therefore  we  also  carried  out  research  on  how  the  morphology  selection  might  be  affected 
by  the  microscopic  dynamics  of  surface  diffusion  and  attachment  of  molecules  on  the 
crystal-vapor  interface  in  general.  We  have  made  good  progress  in  this  research,  and 
believe  the  ideas  we  introduced  there  will  play  an  important  role  in  any  future  research  in 
this  area. 

The  most  striking  feature  of  the  CuPc  molecule  is  its  very  anisotropic  shape:  it  is 
essentially  a  flat  square  plate.  See  Fig.  1  Many  features  of  the  dynamics  and  crystal 
structure  must  be  strongly  influenced  by  this  unusual  shape.  These  shape  considerations 
should  have  more  gener^  implications  for  the  growth  of  organic  thin  films  under  a  variety 
of  conditions.  Similar  ideas  have  proved  very  useful  in  other  contexts.  Thus,  the  structure 
of  liquid  metals  is  well  represented  by  a  fluid  of  hard  spheres  of  appropriate  diameter,  and 
basic  features  of  the  isotropic-to-nematic  phase  transition  in  liquid  crystals  is  controlled  by 
the  anisotropic  shape  of  the  individual  molecules[1.5]. 

As  a  first  attempt  to  incorporate  shape  considerations  into  the  growth  of  CuPc  films, 
we  have  studied  a  simplified  "hard  square"  model.  Growth  is  modeled  by  the  sequential 
attachment  of  planar  square  molecules  whose  only  interactions  are  excluded  volume 
interactions  forbidding  neighboring  molecules  to  overlap.  The  most  recent  results  from  this 


simulation  project  are  summarized  in  Section  B.2. 

Lastly,  we  have  focused  much  effort  on  constructing  computationally  feasible 
simulation  models  of  CuPc  crystal  growth  at  the  molecular  scale,  taking  account  of 
attractive  as  well  as  repulsive  interactions.  To  reach  this  goal  we  have  developed  a 
simplified  intermolecular  potential  and  determined  some  of  its  low  energy  crystal 
structures.  We  have  shown  that  the  simplified  models  are  capable  of  qualitatively 
reproducing  some  basic  features  of  known  CuPc  crystal  structures.  We  have  also  begun  an 
investigation  of  the  energetic  and  kinetics  of  diffusion  of  a  single  molecule  near  a  fixed 
crystalline  substrate.  These  results  are  summarized  in  Section  B.3.  A  detailed  Report  is 
attached  as  Appendix  B. 

B.  Status  of  the  research  effort 

1.  Studying  the  possibility  of  convection  in  Physical  Vapor  Transport  (PVT)  in 
ground-based  and  low  Earth  orbit  ( LEO)  experiments  ( Fisher) 

As  indicated  in  our  last  Report,  we  have  concluded  that  it  is  unlikely  that  significant 
convection  can  arise  as  a  result  of  the  one-dimensional  (vertical)  temperature  gradient  in 
Debe's  experimental  cell.  This  appears  to  rule  out  the  most  obvious  explanation  for  the 
different  behavior  in  LEO  and  the  ground  samples,  and  we  have  not  pursued  this  line  of 
research  much  since  our  last  Report.  We  asked  Professor  Alexander  Chernov,  a  world- 
renowned  authority  on  both  theoretical  and  experimental  aspects  of  crystal  growth,  to 
review  our  conclusions.  He  agreed  with  our  assessment,  but  mentioned  that  radial 
(horizontal)  temperature  gradients  between  the  center  and  sidewalls  of  the  cell  are  also 
present  and  could  conceivably  induce  convection.  However,  a  rough  order  of  magnitude 
estimate  suggests  only  a  very  small  effect.  A  final  Report  detailing  these  conclusions  is 
attached  as  Appendix  A. 

2.  Microscopic  models  of  stochastic  growth  with  contact  interactions  (Einstein) 


Research  continued  on  simulated  crystal  growth  using  the  square  deposition  computer 
program  originally  written  by  L.  Kieffer  Warman,  and  described  in  detail  in  a  previous 
Report.  The  code  produces  an  output  file  consisting  of  the  positions  of  the  centers  of  each 
square  plus  its  orientation.  A  new  square  is  accepted  if  its  normal  makes  an  angle  with  the 
square  it  first  touches  that  is  less  than  some  specified  critical  angle,  with  an  adjustable 
number  (and  distance)  of  lateral  moves  allowed  if  a  trial  exceeds  this  angle.  Histograms  of 
the  orientation  9  of  the  normal  vector  of  each  square,  relative  to  the  vertical  direction,  in  a 
given  pile  were  similar  but  not  identical  to  gaussians,  as  illustrated  in  Fig.  2.  This 
distribution  was  obtained  by  choosing  8  randomly  between  0  and  %.  To  produce  a 
spherically  symmetric  distribution  of  normals,  one  must  choose  a  9  distribution  that  goes 
as  sin  9  .  Histograms  were  made  of  the  piles  produced  with  this  modified  distribution,  with 
quite  different  results:  Instead  of  being  peaked  at  0  degrees,  the  histograms  tended  to  be 
peaked  at  a  value  between  70°  and  80°.  The  spherically  symmetric  distribution  of  square 
normals  is  seen  to  be  a  poor  choice  for  modeling  the  growth  of  CuPc  microcrystals,  which 
have  a  herringbone  structure  with  overlapping  molecules  having  identical  orientations.  The 
spherical  distribution  of  square  normals,  on  the  other  hand,  biases  against  such  behavior, 
because  only  a  relatively  small  number  of  horizontal  squares  are  produced;  the  majority  will 
be  closer  to  the  vertical.  To  model  the  configurations  of  CuPc  more  accurately,  we  turned 
to  a  distribution  that  favored  squares  with  horizontal  (9  =  0)  orientations,  illustrated  in  Fig. 

3. 


Improvements  were  also  made  to  the  deposition  algorithm.  The  version  developed  by 
Warman  and  used  during  the  earliest  work  relied  on  a  discrete  approximation  to  check 
whether  an  incoming  square  made  contact  with  previously  deposited  squares.  Occasionally 
the  code  suggested  that  contact  had  been  made  when  graphical  displays  showed  that  it  had 
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not  As  the  grid  was  made  finer,  the  height  of  a  pile  decreased,  but  it  was  difficult  to  assess 
the  limiting  situation  due  to  limited  memory  and  disk  space.  Hence,  the  program  was 
revised  (with  considerable  labor)  to  include  an  analytic  determination  of  the  point  of  contact 
between  squares. 

From  the  numerical  data  implicit  in  the  piles,  a  major  goal  was  to  characterize  and 
interpret  the  correlations  of  the  orientation  of  squares  as  a  function  of  their  separation,  in 
particular  to  find  the  form  and  characteristic  length  of  the  decay  in  orientation  correlation  of 
pairs  of  squares  as  their  horizontal  separation  increased  with  vertical  separation  fixed,  or 
vice  versa.  Initial  work  did  show  evidence  of  decay  of  orientational  correlations  with 
separation  in  either  direction  and  strong  net  orientational  alignment  for  the  entire  pile.  In  the 
first  stage,  we  generated  scatter  plots  to  illustrate  this  tendency.  Last  summer  we  extended 
this  approach  to  look  at  various  statistical  properties.  The  density  was  checked  as  a  function 
of  height  z  in  the  pile,  as  illustrated  in  Fig.  4.  There  is  a  very  high  density  near  the 
substrate,  where  the  squares  are  nearly  horizontal.  The  density  decays  quickly  to  a  constant 
value  perhaps  half  the  near-substrate  value  (but  dependent  on  the  acceptance  angle)  and 
remains  around  this  value  until  quite  near  the  top  of  the  pUe,  when  it  drops  quickly  to  zero. 
We  examined  the  correlations  of  tilt  angle  by  computing  ((0(r)-9(O))2),  where  r  is  the 
horizontal  distance  from  a  reference  square  (averaged  over  the  entire  pile).  Plotted  vs.  r,  as 
in  Fig.  5,  this  function  quickly  rises,  then  slowly  decays.  The  height  of  the  pile  increased 
linearly  with  time,  suggesting  a  steady-state  growth  process.  We  spent  considerable  effort 
in  examining  the  roughness  w  of  the  pUes,  where  w  is  the  root-mean-square  fluctuations  of 
the  height  of  the  pile,  averaged  over  r.  A  log-log  plot  of  w  vs.  the  average  height  at  first 
grew  linearly,  then  rather  ab^rupdy  leveled  out  at  a  height  that  increased  with  the  lateral 
dimension  L  of  the  pile,  as  shown  in  Fig.  6.  As  illustrated  in  Fig.  7,  we  also  found  that  the 
L-dependence  of  w  could  be  adequately  described  by  the  a  two-parameter  fit  to  the  form  w^ 
=  Cl  -I-  C2*ln(L),  the  behavior  for  the  Edwards-Willdnson  model[2.1],  which  is  expected  at 
relatively  early  times.  (Later  on,  the  growth  paradigm  is  expected  to  shift  to  Kardar-Parisi- 
Zhang[2.2]  for  this  growth  algorithm.)  We  did  note  that  inconclusive  evidence  that  for 
large  L  the  roughness  seemed  to  grow  faster  than  the  above  equation. 

Due  to  the  premature  termination  of  funding,  we  were  not  able  to  analyze  these 
properties  more  systematically,  nor  to  consider  other  properties  of  concern.  E.g.,  we  could 
not  study  how  the  density  and  volume  of  the  piles  depend  on  the  amount  of  diffusion 
(hopping  after  unsuccessful  adsorption  attempts).  We  will  also  did  not  have  time  to  develop 
simple  models  to  compare  with  the  simulated  data.  We  reemphasize  that  the  gener^ 
problem  of  deposition  of  highly  asymmetric  objects  has  received  almost  no  attention,  in 
stark  contrast  to  the  many  studies  of  deposition  of  spheres[2.3],  and  so  seems  a  particularly 
fertile  area  for  future  study. 

3.  Realistic  intermolecular  potential,  crystal  structures,  and  dynamics  of  growth 
(Weeks) 

As  discussed  in  our  previous  reports,  an  empirical  potential  for  the  CuPc  molecule  is 
available  from  the  standard  CHARMM  molecular  modeling  package,  involving  interatomic 
potentials  for  all  the  57  atoms  in  the  CuPc  molecule  [3.1].  Atom-atom  model  potentials  of 
this  kind  have  proved  capable  of  fitting  many  properties  of  large  organic  and  inorganic 
molecules,  and  programs  such  as  CHARMM  are  routinely  used  to  assess  the  feasibility  of 
new  compounds  in  drug  design.  The  CHARMM  potentii  for  CuPc  has  been  shown  to  fit 
one  particular  crystal  structure  of  the  material  to  very  good  accuracy.  However,  this 
potential  requires  a  total  of  57  x  57  =  3249  separate  interaction  terms  to  represent  the  total 
energy  of  just  one  pair  of  CuPc  molecules,  making  large-scale  simulations  impossible  due 
to  the  computation^  load. 

The  first  question  we  addressed  was  whether  it  was  possible  to  develop  a  simplified 
intermolecular  potential  which  preserves  essential  features  of  the  full  potential  while 
requiring  much  less  computation.  The  idea  was  to  represent  groups  of  atoms  in  the  original 


molecule  by  interaction  sites.  We  then  tried  to  choose  empirical  potentials  between 
interaction  sites  in  different  molecules  to  reproduce  as  closely  as  possible  basic  features  of 
the  full  CuPc  intermolecular  potential  as  given  by  the  CHARMM  model.  Particular  attention 
was  paid  to  the  low  energy  "sliding"  configurations  where  the  planar  molecules  are  parallel 
and  the  centers  move  relative  to  one  another.  These  configurations  should  be  very 
important  for  surface  diffusion. 

Clearly  there  is  a  tradeoff  between  an  accurate  representation  of  the  full  potential 
surfaces  and  a  simple  representation  using  only  a  few  interaction  sites.  We  do  not  know 
how  to  quantify  this  tradeoff  and  can  only  use  our  best  judgment  as  to  when  an  acceptable 
compromise  has  been  achieved.  In  our  final  version,  we  introduced  a  model  with  13 
interaction  sites,  not  all  of  which  interact  with  each  other.  With  this  potential,  the  total 
number  of  terms  needed  to  calculate  the  interaction  of  two  molecules  is  brought  down  to 
89.  The  computer  time  in  typical  applications  was  reduced  by  nearly  a  factor  of  100. 
Memory  requirements  are  also  reduced  by  more  than  3/4.  A  picture  of  the  "model" 
molecule,  superimposed  over  the  CuPc  molecule,  is  shown  in  Fig.  8. 

Minimum  energy  crystal  structures  for  this  model  potential  were  calculated  using  a 
multi-dimensional  minimization  code  [3.2].  Multiple  minimizations  from  different  starting 
configurations  yielded  several  different  crystal  structures  that  aU  had  energies  per  unit  cell 
within  a  few  percent  of  the  value  for  the  full  potential.  CuPc  is  known  to  have  several 
crystalline  polymorphs.  Unfortunately  there  is  almost  no  information  in  the  literature  about 
the  details  of  these  different  crystal  structures.  Thus  we  were  not  able  to  determine  whether 
there  is  detailed  correspondence  between  them  and  the  minimum-energy  structures  we  have 
found  using  our  model  potential.  The  main  point,  on  which  both  our  models  and  the 
experiments  agree,  is  that  CuPc  has  a  number  of  different  locally  stable,  low  energy  crystal 
structures,  including  those  with  a  herringbone  structure. 

We  could  easily  tune  the  parameters  of  our  model  further  to  reproduce  even  more 
features  of  any  particular  experimental  crystal  structure.  However,  our  intention  is  not  to 
accurately  reproduce  one  particular  crystal  structure  but  to  get  an  overall  qualitatively 
accurate  potential  describing  interactions  of  very  anisotropic  plate-like  molecules  in  the 
variety  of  configurations  encountered  during  the  course  of  surface  diffusion  and  crystal 
growth.  For  this  purpose  we  believe  our  present  potential  should  prove  adequate.  This 
work  has  been  submitted  for  publication,  and  is  presently  under  review  at  the  Journal  of 
Crystal  Growth.  A  preprint  of  this  work  is  attached  as  Appendix  C. 

Work  has  begun  on  a  second  project,  and  we  hope  to  continue  this  further  on  our 
own  and  eventually  produce  a  second  publication.  The  idea  is  to  determine  the  energy 
barriers  and  kinetics  associated  with  diffusion  of  a  single  CuPc  molecule  on  the  cleaved 
surface  of  the  CuPc  crystal.  Various  techniques  have  been  described  in  the  literature  [3.3] 
for  mapping  out  the  energy  landscape  for  diffusion  on  a  crystal  surface.  In  our  case  we  do 
not  expect  much  relaxation  of  the  crystal  surface  to  occur  during  diffusion;  thus  we  chose 
to  use  a  simplified  technique  in  which  surface  relaxation  is  neglected.  However  the  extra 
degrees  of  freedom  due  to  the  non-spherical  shape  of  the  CuPc  molecule  still  make  this  a 
challenging  calculation. 

We  then  plan  to  use  transition-state  theory  [3.3]  to  estimate  the  kinetics  of  diffusion 
on  this  surface.  We  believe  that  surface  diffusion  may  proceed  in  a  very  different  way  for 
such  large  anisotropic  molecules  than  the  usual  hopping  mode  adopted  by  atoms,  since 
rotational  degrees  of  freedom  must  now  play  a  significant  role.  Debe  has  rightly 
emphasized  the  crucial  role  surface  diffusion  must  play  in  CuPc  crystal  growth,  and  further 
research  along  these  lines  is  definitely  called  for. 


C.  Publications 


A  paper  describing  the  construction  of  our  simplified  model  potential  for  the  CuPc 
molecule  and  the  crystal  structures  it  predicts  has  been  submitted  to  the  Journal  of  Crystal 
Growth  and  is  presently  under  review.  With  some  further  work,  we  hope  to  write  a  second 
paper  on  surface  diffusion  with  anisotropic  molecules.  We  intend  to  write  a  short  paper 
reporting  the  summer  work  on  hard-square  deposition,  but  are  dependent  on  participation 
by  the  student.  The  Final  Report  in  Appendix  A  describing  our  conclusions  with  regard  to 
convection  in  the  cell  has  been  written  in  a  very  detailed  form,  to  document  thoroughly 
what  we  have  done  and  to  aid  any  future  workers  on  this  project.  We  plan  to  publish  a 
short  note  giving  our  main  conclusions,  probably  in  the  Journal  of  Crystal  Growth. 

D.  Personnel  associated  with  the  research  effort 

•  Senior  Personnel:  Professors  T.  L.  Einstein,  M.  E.  Fisher,  and  J.  D.  Weeks. 

•  Postdoctoral  Associate:  Dr.  Robin  Selinger 

•  Graduate  students:  Sheng-nan  Lai  and  Dajiang  Liu 

•  Summer  Student:  David  R.  Eisner 

•  Consultant:  Professor  Alex  Chernov 

E.  Interactions 

•  Dr.  Robin  Selinger  presented  a  paper  at  the  March  1993  meeting  of  the  American 
Physical  Society  entitled  “Simulation  studies  of  morphology  selection  in  the  growth  of 
crystalline  organic  thin  films.”  She  also  presented  a  seminar  on  this  work  at  Georgetown 
University  in  March,  1994. 

•  Professor  Alex  Chernov  read  our  preliminary  technical  report  about  the  possibility 
of  convection  in  the  experimental  cell.  He  visited  the  University  of  Maryland  to  consult 
with  us  about  its  conclusions  and  to  offer  his  suggestions  for  future  work. 
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Figure  Captions 

Fig.  1.  Qiemical  constituents  of  CuPc.  The  nx)lecule  has  a  nearly  square,  plan^ 
structure  with  a  side  length  of  about  10  A  and  a  thickness  of  about  3.5  A. 

Fig.  2;  Histogram  of  the  orientation  q  (relative  to  the  vertical  direction)  of  300, (X)0  squ^es 
in  a  pile,  obtained  choosing  tiie  three  Euler  angles  with  a  flat  random  distribution. 
The  distributions  were  pea^  at  0=^  and  decayed  to  zero  in  a  qualitatively  (but  not 
quantitatively)  gaussian  marmer.  In  addition,  until  0  reached  the  critical  angle  (here 
15°),  the  curve  was  rather  linear,  as  noted  in  earlier  work.  At  the  critical  angle,  a 
discontinuity  in  the  slope  is  evident 

Fig.  3:  Histogram  of  3(X),000  squares  in  a  pile,  obtained  from  the  rnore  appropriate 
distribution  function  ,  where  0  is  measur^  in  radians.  (This  distribution  has  a 
maximum  at  0  =0  and  a  minimum  at  0  =  7c/2  radians.)  As  in  the  histogram  in 
Fig.  1  produced  using  the  original  flat  distribution  of  q  values,  the  histograms  of 
the  square  piles  produced  with  the  new  distribution  function  were  peaked  at  0  =  0 
and  fell  off  towards  zero  as  0  approached  n/l  radians.  Once  again,  disccmtinuity  in 
slope  is  observed  at  the  critical  angle,  here  15°  as  in  Fig.  2. 

Fig.  4:  Density  vs.  z  (coordinate  in  the  growth  direction)  fw  a  1000x1000  deposition  area. 

Fig.  5:  Orientational  correlations  vs.  horizontal  separation,  averaged  over  the  entire  pile. 

Fig.  6:  Log-log  plot  of  roughness  (w)  vs.  mean  height  of  the  pile,  for  several  values  of  the 
length  of  a  side  of  the  [square]  deposition  area,  showing  the  sharp  crossover  to 
saturated  roughness  with  increasing  height 

Fig.  7:  Plot  of  saturated  mean  square  roughness  (w^)  vs.  the  length  of  a  side  of  the 
[square]  deposition  area,  along  with  a  best-fit  to  the  Edwards- Wilkinson  scaling 
form. 

Fig.  8:  The  simplified  molecular  structure  used  in  our  model  potential  for  CuPc.  This 
should  be  compared  to  Fig.  1. 
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IN  THE  3M-PVTOS  EXPERIMENT 


Sheng-Nan  Lai 


ABSTRACT 

Inspired  by  the  surprising  results  found  by  M.  K.  Debe  et  al.  in  their  3M-PVTOS 
experiments,  this  short  report  is  a  preliminary  study  to  understand  the  claimed  impor¬ 
tant  fluid  dynamics  of  physical  vapor  transport  in  a  closed  cylindrical  ampoule.  The 
purpose  and  design  of  the  3M-PVTOS  experiment  are  briefly  discussed.  The  primary 
features  of  films  grown  in  microgravity  environments  and  ground  control  films  are  pre¬ 
sented.  Based  on  the  PVT  parameters  provided  by  M.  K.  Debe  et  al.,  a  crude  estimate 
of  the  possibility  of  convection  is  made.  Finally,  a  summary  of  two  theoretical  studies 
is  presented. 


1.  INTRODUCTION 


Physical  vapor  transport  (PVT)  is  an  experimental  method  for  growing  crystalline 
solids  and  films  by  gaseous  diffusion.  In  the  PVT  process  a  material  is  sublimed, 
frequently  into  an  inert  buffer  gas,  and  allowed  to  advect  down  a  temperature  gradient 
to  a  growth  interface  where  the  material  recondenses.  Growth  of  single  crystals  by 
PVT  ha-s  been  commonplace  for  years.  In  particular,  vapor  transport  within  a  closed 
ampoule  offers  an  experimental  simplicity  often  used.  The  fluid  dynamics  involved 
in  such  a  process  is  not  necessarily  simple,  however.  Buoyancy-driven  convection 
may  arise  in  groimd-based  PVT  experiments,  and  its  presence  certainly  affects  the 
mass  transport  rate  and  temperature  distribution  which  in  turn  affects  the  quality  of 
crystals.  When  the  buffer  gas  is  introduced,  either  by  thermal  decomposition  of  the 
source  material  or  by  intentionally  introduced  inert  gas,  the  interplay  of  convection  and 
diffusion  complicates  the  situation  even  more  and  makes  the  outcome  of  the  experiment 
harder  to  control.  Experiments  in  a  microgravity  environment  offer  an  opportunity  to 
further  study  the  underlying  mechanisms  that  are  influenced  by  unit  gravity. 

In  1985,  M.  K.  Debe  and  co-workers  grew  organic  thin  films  of  copper  phthalocya- 
nine*  (CuPc)  by  PVT  in  both  low  earth  orbit  (LEO)  on  the  Space  Shuttle  Orbiter 
and  in  the  laboratory  as  ground  controls.  Follow  M.  K.  Debe,  we  name  the  experi¬ 
ment  PVTOS  hereafter,  for  PVT  of  organic  solids.  After  a  complete  characterization, 
M.  K.  Debe  et  ai,  in  a  series  of  three  papersj^”^'  reported  that  CuPc  thin  films  grown  in 
LEO  environments,  when  compared  to  those  of  the  ground  control  films,  are  optically 
more  homogeneous  and  smoother,  are  significantly  denser  in  their  central  portions, 
are  more  highly  imiaxiaUy  oriented,  and  contain,  essentially,  a  new  crystalline  poly¬ 
morph  different  from  those  seen  in  ground  control  films.  In  order  to  give  a  satisfactory 
explanation  to  their  surprising  results,  M.  K.  Debe  et  al.  made  a  further  analysis 
of  the  gas-phase  evolution  and  convective  heat  transfer  of  their  PVT  experiments 
and  concluded  that,  among  other  PVT  parameters,  the  presence  of  convection  in  the 
ground-based  PVT  experiments  and  in  the  edge  regions  of  the  LEO  experiments  may 

*  See  Appendix  A  and  Fig.  A.l  for  the  chemical  structure  of  the  CuPc  molecule. 
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be  the  primary  feature  causing  the  small  variation  among  the  ground  control  films  and 
their  larger  distinctions  firom  the  LEO  films/ 

Inspired  by  their  observations,  om  aim  in  this  study  is  to  understand  the  gas-phase 
behavior  of  the  closed-cell  PVT  process  including  the  possibility  of  natural  convection 
(either  thermal  or  solutal)  and,  if  the  convection  exists,  the  flow  pattern,  the  speed 
of  the  convection  and  the  distribution  of  the  temperature  within  the  PVT  ampoule. 
In  particular,  we  would  like  to  estimate  (i)  the  rate  of  arrival  of  CuPc  molecules,  (ii) 
the  speed  of  the  wind  of  residual  CuPc  molecules  and  buffer  gases  sweeping  across 
the  growth  interface,  (iii)  the  temperature  distribution  on  the  growth  interface.  The 
hope  is  that  through  the  understanding  of  this  gas-phase  behavior  one  can  provide  the 
appropriate  microscopic  boundary  conditions  on  the  gas-solid  interface  as  an  input  for 
a  computer  simulation  of  the  crystal  growth.  By  this  approach  we  wish  to  single  out 
the  most  important  parameters  in  the  closed-cell  PVT  process  and,  hopefully,  simulate 
the  advantages,  if  any,  of  the  microgravity  enviromnent  in  an  earth-bound  laboratory. 

In  this  report,  we  first  describe  some  of  Debe’s  experimental  results  and  his  ex¬ 
planations  of  the  surprising  results  (see  also  Appendix  B.)  By  using  the  results  of 
linear  stability  analysis  of  natural  convection  and  its  experimental  verification,  we 
next  analyze  the  possibility  of  natural  convection  in  Debe’s  PVT  experiments.  Then, 
a  summary  of  two  nimierical  studies  in  direct  support  of  Debe’s  PVT  experiments  is 
presented.  Table  III  and  Appendix  A  present  quantitative  estimates  of  many  of  the 
physical  parameters,  macroscopic  and  microscopic  needed  in  assessing  aspects  of  the 
experiments  or  in  setting  up  simulations.  Finedly,  at  the  conclusion  of  this  report,  we 
comment  briefly  on  future  experiments. 


t  See  the  first  paragraph  of  Sec.  5  on  page  345  of  Ref.  3. 
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2.  DISCUSSION 


The  appaxatus  used  by  M.  K.  Debe  et  al.  is  schematically  shown  in  Fig.  1.  Fig.  2  is 
a  cut-away  view  of  the  same  apparatus.  For  a  more  detailed  description  of  the  design 
and  performance  Refe.  4  and  5  may  be  consulted.  The  PVT  process  occurs  within  the 
ampoule,  which  is  a  1.7  cm  diameter  by  7.5  cm  long  Pyrex  tube.  The  organic  source 
material  is  confined  to  the  hot  end  of  the  tube;  the  substrate  platen,  which  is  a  1.4  cm 
diameter  copper  disk,  is  at  the  cooler  end.  Outside  the  ampoule,  there  is  an  electrical 
heater  wound  with  resistance  wire  used  to  maintain  a  non-linear  axial  temperature 
gradient  between  the  hot  end  and  the  cool  end  of  the  ampoule.  A  heat  pipe  maintains 
the  substrate  temperature  near  70°  C;  a  computer-controlled  heater  kept  the  hot  end 
of  the  ampoule  at  about  400°  C  for  4  hours  during  the  PVT  process.  The  ampoule  and 
heater  are  placed  concentrically  into  a  cylindrical  stainless  steel  cell  which  provides  a 
vacuum  envelope  for  the  heater.  At  the  beginning,  the  ampoules  are  fiUed  with  550 
mTorr  He  or  400  mTorr  Xe  and  the  cells,  the  region  outside  the  ampoule  but  inside 
the  stainless-steel  cell  (see  Figs.  1  and  2),  are  evacuated.  After  the  PVT  process, 
the  ampoules  contain,  mainly,  the  original  buffer  gas,  plus  H2,  CO2  and  CO  (of  total 
pressure  about  3  Torr)  which,  except  for  the  original  buffer  gas,  are  outgassed  from  the 
bulk  of  the  glass  and  stainless-steel  components  inside  the  ampoule  and  formed  from 
the  decomposition  of  the  CuPc  source  material.  The  cells’  gas  content  after  the  PVT 
process,  also  due  to  outgassing,  is  primarily  of  H2  with  smaller  amoimts  of  several 
other  common  gases  at  total  pressures  on  the  order  of  one  Torr. 

Fig.  3  and  Fig.  4  show  the  micrographs  of  a  LEO  sample  rind  a  ground  control 
sample,  respectively.  The  groxmd  control  sample  G2,  which  by  virtue  of  the  PVT 
parameters  used,  has  been  considered  by  the  authors  the  best  ground  control:  G2  has 
almost  the  same  PVT  parameters  as  LEOl  has.  They  have  the  same  initial  buffer 
gas,  Xe  at  400  mTorr,  the  same  final  cell  pressure,  0.4  Torr  H2,  and  both  of  them 
have  their  substrate  platens  coated  with  H2PC  seed  films.  The  main  differences  are  the 
temperature  of  the  substrate  platens,  68°  for  LEOl,  77°  for  G2,  and,  of  course,  the 
presence  or  absence  of  gravity,  the  LEOl  sample  was  produced  in  the  microgravity  en- 
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vironinent  while  G2  was  made  in  the  ground-based  laboratory  with  the  ampoule’s  hot 
end  up.  However,  when  one  compares  these  two  sets  of  micrographs,  the  distinctions 
between  the  LEO  sample  and  G2  sample  are  obvious.  Firstly,  the  central  portions  of 
the  LEO  film  is  denser  than  the  central  portion  of  the  ground  control  film  and  their 
microstructures  are  also  different  in  that  region.  Secondly,  the  edge  region  of  the  LEO 
film  is  similar  to  the  edge  region  of  the  ground  control  film. 

Being  so  similar  in  their  PVT  parameters  used,  the  dramatic  difference  between  the 
film  grown  in  ground-based  experiment  and  that  of  the  LEO  is  striking.  One  possible 
cause  of  the  difference,  as  strongly  suggested  by  M.  K.  Debe  et  ai,  is  buoyancy-driven 
convection.  According  to  the  authors,  the  radial  variation  of  the  morphology  of  the 
LEO  films  is  the  consequence  of  convection  winch  is  minimized  over  the  centers  of 
the  LEO  films  while  remaining  significant  at  their  edges  where  very  strong  thermal 
gradients  persist  between  the  platen  edges  and  the  ampoule  walls:  see  Figs.  1  and  2. 
In  the  LEO  film  centers,  where  the  convection  is  minimized,  the  film  growth  mecha¬ 
nism  remains  imchanged  as  the  ampoule  pressure  increases  and  the  film  morphology 
continues  to  exhibit  dense  domains,  resembling  the  base  layer  microstructure.  On  the 
edges  of  the  LEO  films,  and  in  all  regions  of  the  ground  control  films,  the  onset  of 
convection  is  believed  to  cause  a  significant  change  in  the  film  growth  mechanism, 
resulting  in  the  observed  discrete  whiskered  microstructure.  Debe’s  main  arguments 
identifying  convection  as  the  causative  factor  are  presented  in  Appendix  B  in  his  own 
words. 

Now  buoyancy-driven  convection  is  a  result  of  imbalance  of  forces.  The  forces  can 
be  analyzed  in  the  following  way.  Consider  a  bubble  of  warm  fluid,  of  radius  R,  which 
moves  upward  with  a  constant  speed  u  in  a  fluid  where  a  constant  Unear  temperatme 
gradient  is  maintained.  The  bubble  coming  from  the  hotter  region  is  now  in  a  cooler 
environment  with  fluid  density  greater  than  the  fluid  density  in  the  bubble.  The 
difference  in  density  gives  rise  to  the  Archimedean  buoyancy  force,  Fjj  g  \  5p  \  R^. 
Defining  —dp/pdT  =  a  to  be  the  thermal  expansivity  of  the  fluid,  the  buoyancy  force 

*  See  the  first  paragraph  of  Sec.  5  starting  at  page  345  of  Ref.  3. 
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can  be  written  in  terms  of  the  temperature  difference  6T  as  Fj,  gap  \  ST  |  R^. 
The  relaxation  time  of  the  bubble’s  temperature,  due  to  heat  diffusion,  is  r  ~  i2^/ «, 
where  k  is  the  thermal  diffusivity  of  the  fluid.  This  means  that  at  a  given  time  t  the 
temperature  within  the  bubble  is  about  that  of  its  surroundings  at  an  earlier  instant 
of  time,  t  —  r;  so  that  at  time  t  the  temperature  difference  ST  between  the  bubble 
and  its  neighborhood  is  ST  ~|  VT  |  vt  ~|  VT  |  R^v/k,  which  yields  the  Archimedean 
buoyancy  force  Fb  ~  gap  \  VT  |  R^v/k.  The  buoyancy  force  is  opposed  by  viscous 
drag,  Fdt  which,  in  a  first  approximation,  is  of  order  T/iit;,  where  r/  is  the  dynamic 
viscosity  of  the  fluid.  Thus,  the  motionless  steady  state  becomes  unstable  when  the 
buoyancy  force  just  overcomes  the  viscous  drag.  The  ratio  of  buoyancy  force  to  viscous 
drag  is  a  dimensionless  number  called  the  Rayleigh  number 

Ra  =  5a  I  Vr  I  (2-1) 

where 

V  =  -q/p,  (2.2) 

is  the  kinematic  viscosity.  The  Rayleigh  number  increases  with  the  size  of  the  bubble; 
the  instability  should  first  appear  for  bubbles  of  maximum  size,  say  the  radius  of  the 
ampoule  a,  compatible  with  the  boundaries  of  the  container.  With  this  rough  argument 
we  cannot,  however,  estimate  quantitatively  the  value  of  the  critical  Rayleigh  number 
for  the  onset  of  convection.^  We  merely  single  out  which  dimensionless  combination  of 
parameters  is  most  relevant  in  the  stability  analysis. 

An  exact  mathematical  description  of  convection  is  very  difficult  and  an  exact 
solution  is  still  lacking  even  for  the  simplest  system  undergoing  vigorous  convective 
motion.  Various  approximate  methods  do  exist,  which  lead  to  an  analysis  of  the 

t  If  we  use  more  accurate  expressions  for  the  Archimedean  force,  F},,  the  relaxation  time,  r,  cind  the 
viscous  drag,  Fd-  Fb  ircg  |  5p  \  fJ®/3,  r  ~  F^jAn,  Fd  ^  QntjvR^  which  is  the  Stokes’  law,  the 
condition  for  the  onset  of  convection,  FbIFd  >  1,  leads  to  Ra  =  ga  |  VT  |  R*/uk  >  18.  However, 
we  should  not  take  the  above  heuristic  argument  too  seriously.  An  instability  criterion  must  be 
derived  from  a  global  consideration  of  the  fluid,  i.e.  ,  the  overall  internal  energy  released  by  the 
buoyancy  force  must  overcome  the  total  energy  dissipated  by  viscosity  and  thermal  conductivity. 


6 


stability  of  various  inodes  of  motion  in  the  fluid  and  hence  to  predictions  of  which 
modes  are  the  most  likely  to  be  observed.  The  results  £ire  approximate,  but  in  some 
instances  the  approximation  is  a  close  one.  For  a  vertical  cylindrical  container  with 
walls  of  infinite  thermal  conductivity  and  an  aspect  ratio  7  =  X/a  =  6,  convection  was 
found  experimentally  to  set  in  near  Ra*  =  220  ±  12f*'  A  theoretical  study,  assuming 
axisymmetry,  gave  Ra*  =  230f’^  Perfectly  insulating  lateral  walls  reduce  the  threshold 
to  Ra*  =  67.4,  according  to  theory!*' 

Based  on  the  data  provided  by  M.  K.  Debe  et  a/.,'*'  we  estimate  various  ther¬ 
mophysical  parameters  in  Appendix  A.  The  Rayleigh  number  we  obtain  is  Ra  = 
(1.1  ±  0.6)  X  10“^,  far  below  the  threshold  apparently  needed  to  cause  a  buoyant  con¬ 
vection.  Even  if  the  sidewall  temperature  is  approximately  isothermal  from  the  hot 
end  to  a  distance  (say,  0.1  of  the  ampoule’s  length)  close  to  the  cool  end,  near  which 
the  temperature  drops  rapidly  and  the  convection  is  confined  to  that  part  of  the  am¬ 
poule,  the  estimated  Rayleigh  number  would  rise  only  to  0.11  ±  0.06,  stiU  well  below 
the  threshold. 

Another  consideration  should  be  mentioned  here,  however.  For  a  two-component 
fluid  with  molecular  weights  of  the  somrce  material  or  solute  significantly  larger  than 
that  of  the  buffer  gas  or  solvent,  which  is  the  case  of  Debe’s  experiments,  density 
gradients  due  to  the  varying  concentrations  of  the  two  components  may  result  in  con¬ 
vection.  For  example,  with  the  hot  end  up,  the  upward  movement  of  a  fluid  bubble 
carries  it  to  a  region  of  hotter  and  “saltier”  {i.e.  ,  higher  solute  concentration)  en¬ 
vironment.  Suppose  the  temperature  in  the  bubble  is  relatively  quick  to  equilibrate 
with  its  surroundings  while  the  solute  concentration  x  is  significantly  slower  to  do  so. 
The  bubble  then  finds  itself  at  about  the  same  temperature,  but  with  less  solute  and 
therefore  less  density,  than  its  neighborhood,  and  so  the  upward  movement  is  rein¬ 
forced  by  buoyancy  furthering  a  convection.  The  relative  importance  of  these  effects 
is  measured  by  another  Rayleigh  number,  the  so-C£illed  solutal  Rayleigh  number 

Rs  =  I  Vx  I  u^/vDab,  (2.3) 

where  =  dp/ pdx-  Clearly  plays  the  role  of  a  while  Dab-,  binary  diffusion 
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coefficient,  plays  the  role  of  /c  in  the  thermal  convection  analysis.  While  the  theo¬ 
retical  analysis  of  solutal  convection  in  a  closed  cylindrical  ampoule  is  not  available, 

^9] 

experimentally  solutal  convection  was  found  to  set  in  at  Rs  =  180  ±  14. 

By  taking  0  =  1,  which  is  a  good  approximation  for  a  low  concentration  binary 
fluid,  we  find  the  solutal  Rayleigh  number  in  Appendix  A  to  be  Rs  =  (2.2±0.8)  x  lO"**. 
With  such  a  small  solutal  Rayleigh  nmnber,  the  solutal  effect  would  seem  to  play  an 
insignificant  role  in  Debe’s  PVT  experiments. 

Dr.  Chernov*  pointed  out  that  at  the  comer  formed  by  two  walls  of  different 
temperatures  (the  comer  boimd  by  the  copper  platen  and  the  ampoule’s  side-wall  of 
Debe’s  PVT  experiment  seems  to  be  of  this  sort),  thermal  stress  may  lead  to  convection 
even  though  the  Rayleigh  numbers  of  the  system  are  small.  Unfortunately,  we  do  not 
know  how  to  estimate  the  importance  of  this  effect. 

Thermal  creep  is  another  possible  source  that  drives  convection.  Along  a  non- 
isothermal  waU  there  is  a  thin  layer  of  fluid,  called  the  Knudsen  layer,  of  thickness  a 
few  mean  free  paths  that  moves  from  the  cold  to  the  hot  region  with  the  (z-component) 
velocity  which  can  be  expressed  as:^‘°' 


^  v  Ail  /rt  j\ 

V  ci  0.8i/— T — .  (2.4) 

az 

In  an  ampoule  used  to  grow  CuPc,  a  typical  set  of  parameters  is  z/  35  cm^/sec,  with 
a  temperature  difference  of  about  673  —  343  K  =  330  K  maintained  over  a  side-wall 
length  L  of  about  7.5  cm.  Under  these  conditions  the  thermal  creep  velocity  is  nearly 
2.4  cm/sec.  It  is  claimed  that  this  cold-to-hot  flow  along  the  side  walls  will  correspond 
to  a  hot-to-cold  core  flow  of  comparable  velocity,  and  will  probably  not  be  negligible 
compared  to  typical  buoyancy-induced  velocities  (see  Ref.  11,  p.l762  first  column.) 
Unfortunately,  we  do  not  know  how  to  incorporate  this  effect  into  our  analysis. 

We  conclude  this  section  by  stating  that  it  is  not  convincing  to  say  that  there 
is  a  buoyancy-driven  convection  in  Debe’s  PVT  experiments.  Furthermore,  since  the 

*  A.  A.  Chernov,  private  communication. 


8 


gravity  is  reduced  dramatically  to  between  10“^  and  10“®  times  the  earth  value  of 
g  =  9.8  m/s^  in  the  microgravity  environment,  buoyant  force  is  dramatically  reduced, 
so  is  the  chance  of  buoyancy-driven  convection.  Although  one  cannot  exclude  other 
possible  somces  of  convection,  it  is  beyond  the  powers  of  the  author  to  establish  their 
importance  or  unimportance. 

3.  REVIEW  OF  PREVIOUS  NUMERICAL  STUDIES 

Many  theoretical  attempts  have  been  made  to  understand  the  fluid  dynamics  in 
the  closed-cell  PVT  process of  which  Ref.  15  by  T.  L.  Miller  and  Ref.  16  by 
D.  E.  Rosner  and  D.  E.  Keyes  are  numerical  studies  in  direct  support  of  Debe’s  PVT 
experiments.  In  this  section  I  shall  summarize  some  of  their  results. 

Using  the  notation  set  out  in  Table  I  (see  page  29),  Miller  employed  the  following 
governing  equations  to  explore  the  diversity  of  possible  fluid  dynamics  with  the  PVT 
ampoule: 

Conservation  of  radial  momentum: 

=  +  V  (V^t;  -  -)  ,  (3.1) 

dt  Poor  \  r/ 

Conservation  of  axial  momentum: 

^  =  -y .  Vw  -  — 1^  +  (3.2) 

dt  podz  \po  J 

Conservation  of  internal  energy: 

^  =  -y .  Vr  -1-  (3.3) 

dt 

Conservation  of  mass  for  component  A: 

^  =  -V-^X  +  Dab'7\,  (3.4) 

where  incompressibility:  V  •  y  =  0,  and  axisymmetry  are  assumed  and  the  condition 
of  conservation  of  total  mass:  dp/dt  -  -pV  •  V  =  0,  has  already  been  embedded  in 
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the  above  equations.  The  constitutive  relation  is  taken  as 


/?  =  /£>0  [1  -  Q!(T  —  To)  +  Pix  —  Xo)]  •  (3-5) 


The  boundary  conditions  on  temperature  are  assumed  constant  on  each  of  the  end 
walls,  while  sidewalls  are  either  adiabatic  (t.e.  insulating)  or  have  a  fixed  temperature 
profile.  Non-slip  conditions  on  the  botmdaries  are  assumed  for  the  tangential  velocity 
components.  Impervious  sidewalls  are  assumed,  and  the  end  walls  have  the  boimdary 
condition  on  the  axial  velocity  component: 


Dab  dx 

w  =  - - 

l-x9z 


The  reason  that  this  particular  form  is  used  instead  of 


(3.6) 


j  = 


-pDab 


dz' 


(3.7) 


is  that  the  former  is  the  mass  average  velocity  relative  to  a  fixed  coordinate  while  the 
latter  is  the  mass  fiux  of  component  A  relative  to  the  mass  average  velocity  of  the 
whole  fimd.  For  a  full  accoimt.  Ref.  17  may  be  consulted. 

There  exists  an  one-dimensional  steady  solution  to  these  equations.  If  no  variations 
in  the  radial  direction  are  allowed  (and  the  side- walls  are  ignored),  the  axial  velocity, 
temperature,  and  solute  profiles  are  given  by:^‘*' 


W  = 
T{z)  = 
X{z)  = 


The  axial  velocity  amplitude,  W,  is  used  as  the  velocity  scale  in  defining  many  di¬ 
mensionless  parameters  (see  Table  II).  The  one-dimensional  solution  is  also  used  in 
computations  as  the  initial  state  with  which  the  ampoule  system  is  marched  forward 
numerically  in  time  until  a  steady  state  solution  is  reached. 
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Before  presenting  Miller’s  computational  results,  a  few  points  regarding  the  PVT 
parameters  used  by  Miller  must  be  mentioned: 

(1)  In  general,  the  substrate  acts  as  a  sink  for  the  source  material.  This  will  drive 
a  TnA.«ig  flow,  called  the  Stefan  flow,  toward  the  substrate  platen.  The  Stefan 
flow  is  significant  only  when  the  partial  pressure  of  the  source  material  is  larger 
or,  at  least,  comparable  to  the  partial  pressure  of  the  buffer  gas.  Miller  did  his 
simulations  with  the  partial  pressure  of  the  source  material  larger  than  that  of 
the  buffer  gas.  However,  the  averaged  partial  pressure  of  CuPc  vapor  in  Debe’s 
experiments  is  about  10“^  Torr,  which  is  much  smaller  than  the  buffer  gas  partial 
pressure,  4  —  8  Torr. 

(2)  Given  a  small  Rayleigh  number  of  order  1/100  in  the  Debe’s  PVT  experiments, 
absolutely  no  buoyant  effect  can  be  anticipated.  Miller  used  a  larger  Rayleigh 
number  of  order  of  1000  in  his  numerical  simulation  on  the  thermal  convection 
for  the  pm^jose  of  illustration. 

(3)  In  his  simulation  of  the  thermosolutal  convection.  Miller  assumed  that  the  diffu- 
sivity  of  the  solute  is  much  larger  than  that  of  heat.  On  the  contrary,  the  actual 
thermal  diffusivity,  «  =  113  cm^/s  is  much  larger,  than  the  binary  diffusion  coef¬ 
ficient,  Dab  =  12.4  cm^/s  in  Debe’s  experiments  (see  Table  III.) 

Fig.  5  shows  the  computed  effects  of  Peclet  number,  Pe  =  LW/Dab,  on  the  mass 
profile  on  the  crystal  interface  under  zero  gravity.  The  Peclet  number  is  a  measure  of 
the  relative  strength  of  advective  and  diffusive  fiuxes.  The  mass  fiux  is  normalized  to 
the  standard  one-dimensional  flxix  used  as  starting  condition.  The  effect  is  monotonic 
in  Pe,  with  increasing  radial  variation  as  the  Peclet  number  increases.  However,  there 
is  very  little  effect  on  the  radially  integrated  mass  flux. 

Miller  also  did  simulations  on  single-component  fluids  under  unit  gravity.  He 
considered  the  case  in  which  the  ampoule  is  vertical,  with  the  hot  end  down,  and 
the  density  difference  due  to  differing  molecular  weights  of  the  material  is  negligible, 
that  is,  no  solutal  effects  on  the  buoyancy  are  considered.  Even  though  Miller  noticed 
that  “the  actual  physical  parameters  as  estimated  by  Dr.  Mark  Debe  of  3M  result 
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APPENDIX  A  Determination  of  Physical 
Parameters  in  Debe’s  Experiments 

In  the  hope  of  getting  a  feeling  for  Debe’s  3M  PVTOS  experiments  on  the  one  hand 
and  providing  basic  information  for  computer  simulations  on  the  other,  we  estimate 
some  useful  thermophysical  parameters  in  this  appendix.  Although  Debe  has  given 
some  parameters  in  one  of  his  papers,*  no  tmcertainties  are  quoted.  This  appendix  is 
self-contained  as  regards  references. 

Some  information  about  CuPc  molecules 

(A.l)  Molecular  structure  of  CuPc  molecules 

The  formula  for  CuPc  molecules  is  CuPc  =  CUC32N8H16.  As  regards  the  molecular 
structure,  please  see  Fig.  Al. 

(A.2)  Molftcnlar  dimensions  of  CuPc  molecules 

The  copper  phthalocyanine  (CuPc)  is  a  rectangular-like  molecule  of  dimensions 
a  X  a  X  6  with  a  c-  lOA  and  6  ca  SA,  where  6  is  estimated  from  the  interaction  energy 
(see  Fig.  A2)  between  two  CuPc  molecules  with  one  on  top  of  the  other. 

(A.3)  The  crystalline  structure  of  CuPc  crystals 

There  are  at  least  seven  distinguishable  crystalline  polymorphs  of  CuPc.  Figs.  A3 
and  A4  present  the  most  frequently  encounted  crystalline  polymorphs  of  CuPc: 
a— CuPc  and  jS—CvPc. 

B.  Estimates  of  thermophysical  parameters 

The  calculations  are  based  on  the  following  gas  composition: 


Species 

Vi 

Mi(g/mole) 

<Ti(A) 

et/ksOi) 

fci(mW/cm-K) 

H2 

■AVM 

2 

2.827 

59.7 

2.65 

CO 

28 

3.69 

91.7 

0.35 

CO2 

44 

3.941 

195.2 

0.322 

Xe 

12^ 

131.3 

4.047 

231 

0.087 

*  See  M.  K.  Debe  et  al.,  J.  Vac.  Sd.  Techno].,  A8  (1990) 
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where  <t  and  e/fc^  are  taken  from  Appendix  B  of  R.  C.  Reid  et  ai,  The  Properties  of 
Gases  and  Liquids,  fourth  edition  (McGraw-Hill,  New  York  1987),  p.733-734;  ki,  the 
thermal  conductivity,  is  taken  from  M.  K.  Debe  et  oL,  J.  Vac.  Sci.  TecbnoL,  A8  (1990), 
p.58.  The  total  gas  pressure  is  p  =  3.2  Torr  at  T  =  300K  and  the  mean  molecular 
weight  is  M  =  ViMi  =  37.2  g/mole. 

(B.l)  Mean  gas  density,  p 


p  =  nM,  M  =  37.2  g/mole  ~  6.2  x  10  g, 
_N  p 
"'~V~kBT‘ 

At  r  =  300  "K 


(A.1) 

(A.2) 


p  =  3.2  Torr  =  (3.2  Torr)(1.33  x  10^  N/m^  •  Torr), 
(3.2  Torr) (1.33  x  10^  N/m^  •  Torr) 

”  ~  (1.38  X  10-23  J/K)(300  K) 

1.03  X  10^^  l/cm^ 

p  =  nM  ~  (1.03  X  10^^  l/cm3)(6.2  x  lO'^^  g) 

~  6.4  X  10~®  g/cm^. 


~  1.03  X  10“  1/m' 


(A.3) 

(A.4) 

(A.5) 


which  agrees  with  what  Debe  obtained:  6.4  x  10“®  gm/ cm  .  Calculations  of  p  from  the 
partial  pressure  data  given  by  M.  K.  Debe  et  o/.(Figs.  5  and  6  of  J.  Vac.  Sci.  TecbnoL, 
A8  (1990)  p.49-60)  show  that  p  is  between  6  x  10“®  g/cm^  and  1  x  10“®  g/cm^. 
Therefore,  a  reasonable  estimate  for  p  is 


p  =  (8  ±  2)  X  IQ-®  g/cml 


(A.6) 


(B.2)  CuPc  vapor  pressure 


t  See  M.  K.  Debe  et  al.,  J.  Vac.  Sci.  TecbnoL,  A  8,  1990,  p.58,  first  column. 
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The  CuPc  vapor  pressure  can  be  obtained  from  the  empirical  formula 


t 


p  =  760  X  Torr,  {A.7) 

where  A  =  13661  ±  253  and  6  =  14.315  ±  0.367.  At  the  hot  end  where  T  =  400  “C  = 
673  K,  we  have 


p  =  ^760  X  1o(-13661/673+14.315)  ^ 

~  (7.9  X  lO"'*  ±  A)  Torr, 

1/2 

(A.8) 

A  =  (7.9  X  lO"'^  Torr)  [(253)^  (In  10/673)^  +  (0.367)2(ln  10)^ 

~  9.5  X  10"‘‘  Torr, 

(A.9) 

p  =  (7.9  ±  9.5)  X  10“^  Torr. 

(A.10) 

At  the  cold  end  where  T  =  70“C'  =  343iir 


P  =  ^760  X  io(-13661/343+14.315)  ^  ^orr 
a  (2,3  X  10“^^  ±  A)  torr, 

A  =  (2.3  X  10"^  Torr)  [(253)^(ln  10/343)^  +  (0.367)2(lnl0)2]^^^ 

~  4.4  X  10“^^  Torr,  (A.  12) 

p  =  (2.3  ±  4.4)  X  10~^^  Torr.  (A.  13) 

By  using  the  same  formula,  Debe  obtains  the  vapor  pressure  of  CuPc  at  T  =  400°  C 
to  be  3.8  X  10“^  Torr.^ 

(B.3)  Mean  free  path  of  buffer  gas  and  CuPc  molecules. 


f  See  D.  Bonderman  et  ai,  Journal  of  Chemical  and  Engineering  Data,  15  (3)  (1970)  p.399 
§  See  M.  K.  Debe  et  ai.,  J.  Vac.  Sd.  Tecbnol.,  A8,  1990,  p.58,  first  column. 
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Using  the  formula  for  the  mean  free  path,/:  ^ 

ruT 

(A14) 

where  n  is  the  particle  mmiber  density  and  er  is  the  cross  section,  and  the  ideal  gas 

law: 

p  =  nfc^r, 

(A.15) 

we  have 

,  firksT 

V  8  <rp  ‘ 

(A.16) 

Expressing  T  in  K,  p  in  Ton,  <t  in  A,  /  in  cm,  and  substituting  ks  =  1.381  x 
10-23  Joule/“K,  1  Ton  =  1  mmHg  =  1/75.006  N/cm^,  lA  =  lO"®  cm  into  Eq.(A.16), 
we  have 

I  =  (6.5  X  10“^)  —  cm.  (>1*17) 

'  ap 

For  the  buffer  gas  at  T  =  500  K,  p  =  5  Ton,  o-  =  10  A,  the  mean  free  path  between 
buffer  gas  is 

/  ~  6.5  X  10”®  cm.  (A.  18) 

If  there  is  no  buffer  gas,  then  at  T  =  673  K,  p  =  7.9  x  lO”"^  Ton,  the  mean  free  path 
between  CuPc  molecules  is 

I  c:i  11.1  cm,  (A.19) 

where  cr  =  50  A^  is  assumed.  From  eq.(A.18)  and  eq.(A.19),  we  conclude  that  CuPc 
molecules  collide  mostly  with  the  buffer  gas.  If  we  take  p  =  Pbuffer  =  5  Ton,  T  =  500  K, 
and  cr  =  50  A2,  then  the  mean  free  path  of  a  CuPc  molecule  is 

/  =  1.3  X  10"^  cm.  (A.20) 

Debe  does  not  give  estimates  of  the  mean  free  path. 

^  See  Kerson  Huang,  Statistical  Mechanics,  second  edition  (John  Wiley  &  Sons  1987),  p.94,  eq.(5.5) 
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(B.4)  The  rate  of  arrival/bombardment  from  vapor  of  CuPc  on  the  surface 

Assume:  (i)  mass  transport  by  diffusion  only,  (ii)  the  buffer  gas  is  stagnant.  The 
mass  flux  of  CuPc,  j,  satisfles  the  equation: 


j  =  -pDab 


dz’ 


(A.21) 


where  x  is  fractional  weight  concentration  of  CuPc  vapor,  Dab  is  the  binary 
diffusion  coeflicient,  and  p  is  the  mean  density  of  the  fluid.  By  deflnition,  j  is  the  mass 
flux  of  CuPc  relative  to  the  mass  average  velocity  of  the  fluid.  For  a  dilute  S3^tem, 
X  1,  and  since  the  buffer  gas  is  stagnant,  j  can  be  regarded  as  the  mass  flux  of 
CuPc  relative  to  a  fixed  coordinate,  e.g.,  the  ampoule.  If  we  assume  that  p  and  Dab 
are  constant  throughout  the  ampoule,  then 


dz 


—  =  constant, 
pDab 


(A.22) 


where  j  =  const,  is  a  consequence  of  conservation  of  mass.  The  solution  to  eq.(A.22) 
is 


XL- Xo  pal- P Ad 

3  =  pDab — ; —  =  Dab - ? - 1 


(A.23) 


where  XL  and  xo  are  the  values  of  x  at  the  hot  end  and  the  cold  end,  respectively,  and 
PAL  and  paq  are  the  density  of  CuPc  at  the  hot  end  and  cold  end,  respectively. 


Using  the  formula 


_  Mapa 

ksT  ’ 


(A.24) 


where  Ma  =  575.67  g/mole  is  the  mass  of  CuPc  molecules,  pA  is  the  vapor  pressure 
of  CuPc,  we  have 


PAL  ^  1-1  X  10  ®  g/cm^,  (^-25) 

PAo  6.2  X  10”^®  g/cm^,  (A.26) 
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and 

j  ~  1.8  X  10“®  g/cro?  •  s.  {A.27) 

which  compare  well  with  what  Debe  obtained  from  the  growth  rate  of  the  CuPc  films.* 

(B.5)  Kinematic  viscosity,  v 

Using 


(9-3.9) 

(9-4.3) 

(9-5.2) 

(9-5.1) 


„  1.16145 


0.52487 


+ 


2.16178 


(17)0.14874  ^  exp(0.7732T7)  exp(2.438727) 

"  [8(1  -h 


<t>ij  = 


^=E 

t=»i 


ym  I 


/iPoise 


(^4.28) 

(^4.29) 

(^4.30) 

(44.31) 


where  T*  =  T/iei/ks)-  The  above  formulas  are  quoted  from  The  properties  of  Gases 
and  Liquids,  3rd  edition  by  R.  C.  Reid  et  ai.  The  proceeding  equation  numbers  are 
the  equation  numbers  appearing  in  that  book.  At  T  =  (343  K  -I-  673  K)/2  =  508  K, 
the  dynamical  mixture  viscosity  is  found  tohe  p  284  ^Poise.  According  to  the  book 
from  which  the  above  formula  are  quoted,  errors  seldom  exceed  3  to  4  percent  for  the 
estimated  value  of  p.  Therefore,  a  conservative  estimate  of  is 


p  =  284(1  ±  5%)  ^Poise  =  284  ±  14  /xPoise 
=  (284  ±  14)  X  10"®  g/s  •  cm.  (A.32) 


The  kinematic  viscosity: 


u  = 


a 

P 


280  X  10  ®  g/s  •  cm 
8  X  10“®  g/cm^ 


±  A  =  (35  ±  A)  cm^/s. 


*  See  M.  K.  Debe  et  oL,  J.  Vac.  Sci.  TecbnoL,  A  8,  1990,  p.58,  near  eq.(ll). 


(A.33) 
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A  =  (35cmVa)  +(|) 


9cin^/s, 


(^.34) 


1/  =  (35  ±  9)  cm^/s. 


(.4.35) 


this  value  is  diflferent  from  what  Debe  obtained:  130  cm^/s.*  The  formulas  Debe 
quoted  for  calculating  thermophysical  parameters  are  wrong  or  misprinted.^ 


(B.6)  Thermal  diffusivity,  k  =  fc/pcp 


Using* 


(.4.36) 


and  yi,  ki  given  in  table  at  the  beginning  of  this  appendix,  the  thermal  conductivity 
of  the  gas  mixture  is 


k  c-  0.8  mW /cm  •  K  1.9  x  10  ^  cal/cm  •  s  •  K  (.4.37) 

The  thermal  conductivity  we  get  is  slightly  larger  than  what  Debe  obtained:  1.6  x 
10”^  cal/cm  •  s  •  K.^  Error  is  within  3%^  Therefore 


fc  =  1.9  X  10  ^(1  ±  5%)  cal/cm  •  s  •  K 


*  See  M.  K.  Debe  et  aL,  J.  Vac,  Sci.  TedmoL,  A  8,  1990,  p.59,  second  column.  Debe  uses  the 
symbol  rj  for  the  kinematic  viscosity. 

t  In  Ref.  4  p.58  Debe  quoted 

1.16145  0.52487  2.16178 

~  r?  exp(0.7732i;*)  exp(2.4387i;*)  ’ 

[8(1  + Mi/M,)] 

If  the  above  formulas,  instead  of  eqs.  (A.29)  and  (A.30),  are  utilized,  then  u  =  437.64  cm^/s. 

•  See  R.  C.  Reid  et  al,  The  Properties  o!  Gases  and  Liquids,  3rd  edition,  Eq.(10-6.1) 
t  See  J.  Vac.  Sci.  Tedmol,  A8(1990)  p.59,  2nd  column. 

§  See  The  Properties  of  Gases  and  Liquids,  3rd  edition,  p.512. 
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=  (1.9  ±  0.1)  X  10  ^  cal/cm  •  s  •  K. 


(A.38) 

The  averaged  heat  capacity  per  unit  mass  as  given  by  M.  K.  Debe  et  al^  is  Cp  = 
0.218  cal/g  •  K.  For  an  ideal  gas  with  mass  M  =  37.2  g/mole 

~  3  R  ^  5.96  cal/mole  •  K  ^  0.16  cal/g  •  K  (^4.39) 

Therefore,  a  reasonable  estimate  for  Cp  is 

Cp  =  (0.21  ±  0.05)  cal/g  •  K.  (^.40) 

For  the  thermal  difiusivity  we  find 

_  A:  _  1.9  X  10“^  cal/cm  •  s  •  K 

PCp  (8  X  10“®  g/cm^)(0.21  cal/g  •  K) 


ci:  (113  ±  A)  cm^/a 

(>1.41) 

A  /110  /O.OSV]^^^ 

A -(113  cm /s)  g  j  ■^(o.2l) 

c-  39  cm^/s. 

(A.42) 

K  =  (113  ±  39)  cm^/s. 

(A.43) 

(B.7)  Binary  diffusion  coefficient,  Dab 

The  value  for  Dab  estimated  by  Debe,  Dab  =  12.4  cm^/s,  is  subject  to  about  20% 
uncertainty: 

Dab  ^  12.4(1  ±  20%)  cmVs  (12.4  ±  2.5)  cmVs,  (>1.44) 

Dab  =  (12.4  ±  2.5)  cmVs.  (^-45) 

(B.8)  Thermal  Rayleigh  niunber,  Ra  =  ga  |  VT  |  o^/uk 

(980  cm/s=^)  (2  x  IQ-^  1/K)  (320  K)  (0.85  cm)^ 

^  (35  cm^/s)  (113  cm^/s)  (7.5  cm) 

^  See  J.  Vac.  Sci.  Tecbnol.,  A8(1990)  p.58. 
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~  1.1  X  10“^  ±  A, 


.  -  [(¥)■*  -  (5S)’  *  (i)’  <m)'- (S) 


(^.46) 

t1/2 


~  0.6  X  10"^  (.4.47) 


Ra  =  (1.1  ±  0.6)  X  10~^ 

(B.9)  Solutal  Rayleigh  number,  Rs  =  [  Vx  |  a^lvDAB 

By  definition 


X  = 


PA  PA 


PA  +  PB  P 


if  we  assume  the  ps  is  constant,  then 


which  leads  to 


(A.48) 


(.4.49) 


(A.50) 


(A.51) 


pdx  1-x’ 

For  a  dilute  system,  x  1>  we  may  take  /?  =  1,  then  the  solutal  Rayleigh  number  is 


(980  cm/s^)  (1.4  x  IQ-^)  (0.85  cm)^ 
^  (35  cm^/s)  (12.3  cm^/s)  (7.5  cm) 


~  2.2  X  10"^  ±  A, 

~  0.8  X  10"^,  (A.53) 


(>1.52) 


1/2 


Rs  =  (2.2  ±  0.8)  X  10 


-4 


(A.54) 


22 


APPENDIX  B.  Summary  of  Debe’s  Arguments  for  Convection 

We  summarize  Debe’s  arguments  for  convecton  in  this  appendix  as  a  handy  refer¬ 
ence. 

“The  strong  radial  dependence  of  the  LEO  films*  apparent  color,  ellipsometric 
response  parameters  and,  most  importantly,  reflection  spectra  suggests  that  near 
the  edges  of  those  microgravity-grown  films,  the  density  is  less  than  at  their  cen¬ 
ters,  while  the  lower  density  of  the  ground-control-grown  films  varies  much  less 
from  one  side  of  the  film  to  the  other.  The  fact  that  the  optical  properties  of 
the  edges  of  LEO  films  are  the  same  as  those  of  the  unit-gravity-grown  films 
provides  a  self-consistent  argument  that  it  is  the  absence  of  the  convection  at 
the  centers  of  the  LEO  films  which  has  produced  the  enhanced  thin  film  proper¬ 
ties  described  in  this  and  the  following  two  papers  of  this  series,  and  not  some 
other  PVT  parameter  or  experimental  artifact.”  “In  ref.  12  evidence  was  pre¬ 
sented  which  indicated  that,  dioring  the  growth  of  the  groimd  control  films,  there 
was  more  convective  heat  transferred  to  the  substrate  end  of  the  ampoule  when 
the  ampoule  was  processed  with  its  hotter  end  below  its  cooler  end  than  vice 
versa.  Also,  those  data  indicated  that  more  convective  heat  transfer  occurred 
when  the  ampoule  was  processed  with  its  hotter  end  up  in  unity  gravity  than 
when  processed  in  LEO.  These  experimental  observations  are  not  unexpected 
given  the  strong  radial  thermal  gradients  known  to  exist  between  the  edges  of 
the  sample  substrate  platens  and  the  adjacent  inside  surface  of  the  ampoule 
walls.  Without  further  consideration  of  other  transport  mechanism  dififerences 
occurring  between  microgravity  and  unit  gravity,  it  would  appear  as  though  the 
presence  of  convection  has  affected  the  microstructure  of  the  vapor-transport- 
grown  films  on  a  submicron  scale  which  is  microscopic  compared  with  the  size 
of  classical  convective  cells.”  —  From  Ref.  1,  sec.  5,  p.286-287. 

“The  existence  of  the  new  polymorph  is  considered  to  be  a  consequence  of  both 
the  closed-cell  PVT  processing  and  the  epitaxial  growth  onto  an  oriented  or¬ 
ganic  substrate  film,  but  its  exclusivity  and  high  degree  of  crystallinity  in  the 
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microgravity-grown  films  are  attributed  to  the  absence  of  buoyancy-driven  con- 
vecton  which  otherwise  may  disturb  the  delicate  conditions  required  for  epitaxy.” 
—  Prom  Ref.  2,  sec.  7,  p.324. 

“Scanning  electron  micrographs  show  that  the  film  microstructure  in  the  centers 
of  the  space-grown  PVTOS  CuPc  is  not  only  dramatically  different  from  the 
ground  controls,  consistent  with  a  new  crystal  form  as  claimed  here,  but  also  at 
least  a  factor  of  2  denser.  This  enhanced  density  is  seen  to  have  a  strong  radial 
character  over  the  diameter  of  the  films,  explaining  the  optical  effects  described 
in  ref.  1  and  attesting  to  the  influence  of  thermal-gradient-driven  convection  on 
the  microscopic  physical  structxue  of  vapor-transport-grown  films.”  —  Prom 
Ref.  2,  sec.  7,  p.324. 

“A  second  observation  is  that  the  parasitic  CuPc  wall  deposition  just  above  the 
platen  (see  Pig.  1  in  ref.  1)  seems  to  be  related  to  the  radial  variation  in  the 
film  areal  density,  mass  per  imit  area,  of  the  LEO  films,  but  less  often  in  the 
ground  control  films.  The  film  areal  densities  are  higher  in  the  centers  of  the 
LEO  films  than  the  edges,  suggesting  that  the  wall  deposition  has  lowered  the 
volume  density  of  the  CuPc  vapor  immediately  above  the  growth  interface  at  the 
edges  to  below  the  voliune  density  above  the  centers  of  the  LEO  films.  Por  the 
same  Tj  at  the  center  and  edges,  the  supersaturation  ratio  of  the  CuPc  vapor 
is  thus  higher  at  the  center  than  at  the  edges  of  the  LEO  films.  In  the  groimd 
control  films,  the  less  obvious  difference  in  areal  mass  density  between  centers 
and  edges  implies  that  convective  mixing  has  lessened  the  variation  in  volume 
mass  density  of  the  CuPc  vapor  over  the  platen  area,  resulting  in  a  more  radially 
uniform  Sr”  —  Prom  Ref.  3,  sec.  4.4,  p.342. 

“The  substrate  temperature  increase  is  larger  for  higher  ampoule  pressures  and 
hot-end-down  runs,  due  to  increased  heat  transfer  from  the  ampoule’s  wall  and 
hot-end  to  the  substrate.  This  clearly  indicates  convection  is  contributing  to  the 
substrate  temperature  in  the  hot-end-down  orientation,  since  thermal  diffusion 
(Soret)  effects  are  probably  second  order  for  heat  transfer  and  not  orientation 
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sensitive.  Given  the  low  Rayleigh  number  for  this  dilute  system  and  ampoule 
aspect  ratio,  a  one-dimensional  axisymmetric  linear  stability  model  would  not 
predict  the  existence  of  any  free  convection.  It  caimot  be  inferred  however  that 
there  is  convection  in  the  hot-end-up  orientation,  cometimes  referred  to  as  the 
thermally  stable  configuration.  Evidence  that  convection  exists  even  in  the  hot- 
end-up  orientation  can  be  obtained  by  using  the  data  of  Figs.  7-9  to  predict 
what  the  substrate  temperatures  of  the  microgravity  processed  substrates  would 
have  been  had  they  been  processed  in  unit  gravity,  and  then  showing  that  such 
temperatures  are  well  above  the  actually  observed  substrate  temperatures.”  — 
From  Ref.  4,  p.55-56. 
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TABLE  I  List  of  Symbols  Used 


a  Radius  of  ampoule. 

Dab  DiflFusivity  of  components  A  and  B  in  the  mixture. 
g  Gravitational  acceleration. 

k  Thermal  conductivity. 

L  Axial  length  of  ampoule. 

p  Pressure  within  fluid. 

g  Heat  current, 

r  Radial  coordinate. 

To  Temperature  at  z  =  0,  the  cold  end. 

Ti,T2  Other  observed  temperatures:  see  Fig.  1. 

Tl  Temperature  a.t  z  =  L,  the  hot  end. 

V  Local  velocity  vector. 

V,  w  Radial  and  axial  velocity  components,  respectively. 

W  One-dimensional  solution  amplitude  for  axial  velocity. 

z  Axial  coordinate. 

a  Thermal  expansivity. 

Solutal  density  factor. 

7/  Dynamic  viscosity. 

K  Thermal  diffusivity. 

u  Kinematic  viscosity. 

p  Fluid  density. 

Po  Reference  density. 

X  Fractional  weight  concentration  of  component  A  (CuPc). 

XQ  Fractional  weight  concentration  at  z  =  0,  cold  end. 

XL  Fractional  weight  concentration  dX  z  =  L,  hot  end. 
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TABLE  II  Definitions  of  Dimensionless  Parameters 


Nu  =  qL/k 
Pe  =  LWIDab 
Pet  =  LW/k 
Pr  =  v/k 

Ra  =  ga(VT)a^/uK 
Re  =  aW/u 
Rs  =  g0{Vx)a^lJ'DAB 
Sc  =  u/Dab 
'r  =  Lla 


Nuaselt  number 
Solutal  Peclet  number 
Thermal  P&iet  number 
Prandtl  number 
Thermal  Rayleigh  number 
Reynolds  number 
Solutal  Rayleigh  number 
Schmidt  number 
Aspect  ratio 
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TABLE  III  Quantitative  Features  of  Debe’s  Experiments 


Average  mixture  kinematic  viscosity 
Thermal  diffusivity^"^ 

Average  heat  capacity^*' 

Binary  diffusion  coefficient^"^ 

Effective  CuPc  vapor  pressure  at  the  growth  interface^*’ 
Vapor  pressure  of  CuPc  at  400® 

Total  gas  pressure  inside  the  ampoule 

Temperature  difference  between  hot  end  and  cold  end^‘^ 

Total  mean  gas  density 


1/  =  (50  ±  9)  cm^/s 
K  =  (113  ±39)  cmVs 
Cp  =  (0.21  ±  0.05)  cal/g  •  K 
Dab  =  (12.4  ±  2.5)  cm^/s 
p  ~  10“^  Torr 
p  ~  8  X  10“^  Torr 
p  ~  4  ~  8  Torr 
AT  =  335  K 

p  =  (8  ±  2)  X  10"®  g/cm^ 


(a)  See  Appendix  A. 

(b)  See  M.  K.  Debe  et  al.,  J.  Vac.  Sci.  Techno/.,  AS  (1990)  49. 
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FIGURE  CAPTIONS 


Fig.  A1  Molecular  structure  of  CuPc  molecules.  (After  M.  K.  Debe  et  al.,  Ref.  2.)  on 
top  of  the  other.  (Communication  from  Mr.  Da-Jiang  Liu.) 

Fig.  A2  The  interaction  energy  between  two  CuPc  molecules  with  one  on  top  of  the  other. 
(Communication  from  Mr.  Da-Jiang  Liu.) 

Fig.  A3  Crystal  structme  of  a-CuPc.  (Communication  from  Dr.  Robin  Selinger.) 

Fig.  A4  Crystal  structme  of  )3-CuPc.  (After  Dr.  Robin  Selinger.) 

Fig.  1.  Schematic  diagram  of  the  demovmtable  PVTOS  ampoule  within  its  cell  envelope, 
showing  details  of  its  internal  components  and  the  three  temperatures  measured, 
Ti  the  hot  end,  T2  the  substrate,  and  the  cell  surface  midpoint  temperature. 
(After  M.  K.  Debe  et  a/..  Ref.  4.) 

Fig.  2.  Cut-away  view  of  a  PVTOS  cell  showing  the  internal  heater  and  ampoule  assem¬ 
bly  used  for  physical  vapor  transport  deposition  of  organic  films  in  the  Orbiter 
mid-deck.  (After  M.  K.  Debe  et  ai,  Rev.  Sci.  Instrum.,  61  (1990)  865.) 

Fig.  3.  Scanning  electron  micrographs  taken  from  the  center  of  the  microgravity  grown 
film  LE02:  (a)  center,0'’;  (b),  (c)  center,  45°;  (d)  edge,  45°.  (Magnifications, 
30000X.)  (After  M.  K.  Debe  et  al.,  Ref.  3.) 

Fig.  4.  Scanning  electron  micrograph  from  the  ground  control  sample  G2  from  two  view¬ 
ing  angles  of  the  center  and  edge  regions:  (a)  center,  45°;  (b)  center,  45°;  (c) 
edge,  45°;  (d)  edge,  0°.  (Magnifications:  (a),  (b),  (d)  30000x;  (c)  lOOOOx.) 
(After  M.  K.  Debe  et  al.,  Ref.  3.) 

Fig.  5.  Normalized  mass  flux  on  the  crystal  interface  as  a  function  of  radial  distance 
from  the  sidewall  for  varying  Pe  and  for  fixed  7  =  1,  Sc  =  1,  and  Pr  =  l,and 
for  no  gravity,  (a)  Pe  =  5.29,  (b)  Pe  =  4.37,  (c)  Pe  =  2.74,  (b)  Pe  =  2.20,  (b) 
Pe  =  0.94.  (After  T.  L.  Miller,  NASA  Technical  Paper  2620.) 
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Fig.  6.  Nusselt  number  (dimensionless  heat  flux)  as  a  function  of  thermal  Rayleigh  num¬ 
ber  (zero  solutal  Ra)  for  the  cases  of  Pe  =  0  (no  crystal  growth)  and  Pe  =  2.2, 
and  for  aspect  ratio  7  =  1.  (After  T.  L.  Miller,  NASA  Technical  Paper  2620.) 

Fig.  7.  Temperature  and  stream  function  contours  for  the  case  with  the  hot  end  down 
and  buoyancy  effects  of  a  heavy  component  A  included  as  well  as  those  due  to 
temperature.  The  parameters  used  are:  aspect  ratio  7  =  2,  Pr  =  1,  Sc  =  0.1, 
Pe  =  1,  Ra  =  3484,  Rs  =  -1742.  (After  T.  L.  Miller,  NASA  Technical  Paper 
2620.) 

Fig.  8.  As  in  Fig.  7,  but  for  Ra  =  —1742  and  Rs  =  6969  (cold  end  down).  (After 
T.  L.  Miller,  NASA  Technical  Paper  2620.) 

Fig.  9.  Grid,  temperature  contours  T  =  510(10)673  K  and  CuPc  vapor  contours  x  — 

1.24  X  10“®(10~^)3.82  X 10“^  at  standard  operating  conditions  (3.2  Torr).  (After 
D.  E.  Rosner  and  D.  E.  Keyes,  NASA  Contractor  Report  185122.) 

Fig.  10.  Grid,  temperature  contours  T  =  510(10)673  K  and  CuPc  vapor  contours  x  — 

1.24  X  10“®(10~'^)3.82  X  10~^  at  standard  operating  conditions  with  parasitic  heat 
loss.  (After  D.  E.  Rosner  and  D.  E.  Keyes,  NASA  Contractor  Report  185122.) 

Fig.  11.  Grid,  temperature  contours  T  =  510(10)673  K  and  CuPc  vapor  contours  x  = 

1.24  X  10”^(10”^)3.82  X  10~^  at  standard  operating  conditions  with  parasitic 
heat  loss  and  equalibriimi  vapor  boundary  conditions  on  the  side-walls.  (After 
D.  E.  Rosner  and  D.  E.  Keyes,  NASA  Contractor  Report  185122.) 

Fig.  12.  Grid,  stream  function  contours,  vorticity  contours,  temperature  contours  T  = 
510(10)673  K  and  CuPc  vapor  contours  x  =  1-24  x  10~®(10”^)3.82  x  10~^  at 
po  =  3.2  X  10^  Torr(=  4.2  atm).  (After  D.  E.  Rosner  and  D.  E.  Keyes,  NASA 
Contractor  Report  185122.) 
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MODELING  ORGANIC  MOLECULAR  THIN  FILMS:  A  REDUCED 
INTERMOLECULAR  POTENTIAL  FOR  COPPER  PHTHALOCYANINE 


Modelling  Orgzmic  Molecular  Thin  Films:  A  Reduced 
Intermolecular  Potential  for  Copper  Phthalocyanine 

DarJiang  Liu,  Robin  L.  Blumberg  Selinger  and  John  D.  Weeks 

Institute  for  Physical  Science  and  Technology,  University  of  Maryland,  College  Park,  MD  20742 

(August  9,  1994) 

Abstract 

The  intermolecular  potential  energy  between  two  molecules  is  often  repre¬ 
sented  as  a  sum  of  spherically  symmetric  atomic  potentials  between  all  pairs 
of  atoms  on  each  molecule.  While  sudi  potentials  encode  mudi  physical  and 
rhpniira.1  information,  the  number  of  terms  that  must  be  summed  over  goes 
as  the  square  of  the  number  of  atoms  in  a  molecule;  this  ineflBciency  restricts 
the  use  of  these  potentials  in  computer  simulation  and  modeling  of  large 
molecules.  Starting  with  such  an  atomistic  potential  for  Copper  Phthalo¬ 
cyanine  (CuPc),  a  molecule  with  57  atoms,  we  determine  a  simpler  reduced 
intermolecular  potential  consisting  of  a  sum  of  effective  pair  interactions  be¬ 
tween  a  much  smaller  number  of  appropriately  chosen  “interaction  sites”  on 
each  molecule.  Some  of  the  general  issues  and  physical  considerations  that 
arise  when  attempting  this  reduction  are  discussed.  We  arrive  at  an  effec¬ 
tive  potential  involving  13  interaction  sites  in  each  molecule.  This  potential 
reproduces  many  qualitative  features  of  the  full  atomistic  potential  model 
for  CuPc  but  is  much  easier  to  evaluate  numerically,  requiring  only  1%  as 
much  computation  time  as  the  full  atomistic  potential.  Crystal  structures  of 
CuPc  using  both  the  atomistic  and  reduced  potentials  are  determined  and 
compared.  Other  possible  applications  of  these  ideas  are  discussed. 
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I.  BSTTRODUCTION 


The  structure  and  properties  of  crystalline  thin  films  of  Copper  Phthalocyanine  (CuPc) 
have  been  the  subject  of  much  investigation^’^  due  to  this  material’s  potential  applications 
in  both  solar  energy  conversion^  and  gas  sensing^.  When  grown  from  the  vapor  phase, 
these  thin  films  show  a  wide  variation  of  morphology  under  changes  in  growth  conditions®. 
Sensitive  dependence  of  morphology  on  growth  conditions  such  as  substrate  temperature 
has  been  observed  in  thin  films  of  other  organic  materials  as  well®. 

The  microscopic  growth  mechanisms  which  determine  the  structure  and  properties  of 
organic  crystalline  thin  films  are  poorly  understood.  One  might  expect  that  in  the  case  of 
CuPc  the  planar,  “tile-shaped”  character  of  the  individual  molecules  plays  an  important 
role  at  the  microscopic  level.  The  CuPc  molecule  has  a  nearly  square,  flat  structure  with 
a  side  length  of  approximately  10  A  and  a  thickness  of  about  3.5  A  ,  as  shown  in  Fig.  1. 
During  crystal  growth  from  the  vapor  phase,  the  fundamental  processes  of  surface  diffusion 
and  attachment  of  such  anisotropic  molecules  could  differ  significantly  from  those  observed 
in  simpler  materials,  i.e.  metal  on  metal  epitaxial  growth^.  In  order  to  make  progress  in 
understanding  these  microscopic  aspects  of  crystal  growth  of  CuPc,  we  must  first  find  an 
intermolecular  potential  function,  which  is  simply  the  energy  of  interaction  between  two 
molecules  as  a  function  of  their  relative  separation  and  orientation.  Through  its  dependence 
on  relative  orientation,  the  potential  will  implicitly  include  detailed  information  about  the 
molecule’s  planar  shape.  The  potential  must  be  efficient  enough  in  computation  time  to  be 
useful  for  modelling  applications. 

An  atomistic  intermolecular  potential  for  CuPc  consisting  of  a  sum  of  pair  potentials 
for  all  constituent  atoms  was  developed  by  A.  Pohorille®  using  potential  functions  and  pa¬ 
rameters  given  in  the  CHARMM  molecular  modelling  software  package^.  Pohorille  adjusted 
some  of  the  relevant  parameters — the  partial  charges  within  the  CuPc  molecule — in  order  to 
achieve  good  agreement  between  the  potential  and  the  known  crystal  structure  of  /3-CuPc. 
As  discussed  below,  such  a  potential  includes  N  x  N  pairwise  interactions  where  N  is  the 
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number  of  atoms  in  a  molecule;  iV  =  57  for  CuPc.  Thus  the  determination  of  the  interaction 
energy  between  just  two  molecules  requires  the  calculation  and  summation  of  57  x  57  =  3249 
separate  pair  interaction  functions.  Such  a  potential  requires  too  much  computation  time 
to  be  useful  for  modelling  applications  involving  more  than  a  few  molecules.  A  look-up 
table  is  not  a  practical  option  because  of  the  large  number  of  degrees  of  freedom;  storage 
requirements  for  the  tabulated  values  would  be  too  large. 

Of  course,  even  this  atomistic  representation  of  the  interaction  is  only  an  approxima¬ 
tion.  Experience  with  such  model  potentials  suggests  that  they  are  capable  of  providing  an 
accurate  representation  of  many  qualitative  and  large-scale  features  of  the  true  intermolec- 
ular  interactions,  particularly  in  a  case  like  this  where  no  strong  intermolecular  bonds  are 
formed  and  the  main  attractive  interaction  is  due  to  van  der  Waals  forces.  However  it  is 
unlikely  that  the  fine  details  of  the  inherently  quantum  mechanical  interactions  are  correctly 
described  using  a  simple  classical  pairwise  additive  model  potential. 

The  question  we  address  here  is:  can  we  derive  a  “reduced”  intermolecular  potential  for 
CuPc  which  reproduces  these  (presumably  correct)  large-scale  features  of  the  full  atomistic 
intermolecular  potential,  yet  is  computationally  efficient  enough  to  use  for  large-scale  sim¬ 
ulation  and  modelling?  We  have  chosen  an  approach  in  which  the  atomistic  representation 
of  the  CuPc  molecule  is  replaced  by  a  relatively  small  number  of  “interaction  sites” ,  which 
can  be  thought  of  as  representing  groups  of  atoms  in  the  original  molecule.  The  reduced 
intermolecular  potential  is  calculated  as  a  sum  of  pair  potentials  summed  over  all  pairs 
of  interaction  sites.  The  individual  pair  potentials  for  different  species  of  interaction  sites 
are  parameterized  in  the  Lennard-Jones  6-12  form,  with  parameters  fitted  to  reproduce  as 
closely  as  possible  important  features  of  the  full  atomistic  intermolecular  potential.  In  this 
way,  we  were  able  to  determine  a  reduced  potential  which  reproduces  the  gross  features  of 
the  atomistic  potential,  yet  requires  only  1%  as  much  computation  time. 

An  alternative  approach  is  to  represent  the  molecular  interaction  as  a  single  site-site 
potential  function  which  depends  explicitly  on  the  center-to-center  separation  between 
molecules  and  their  relative  orientation.  This  approach  was  first  proposed  by  Corner^®, 
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Gay  and  Berne^^  showed  that  for  a  rigid  linear  molecule  consisting  of  four  Lennard-Jones 
particles,  the  16  atom-atom  pair  interactions  could  be  well  represented  by  a  single  site-site 
potential  defined  as  the  overlap  between  ellipsoidal  Gaussian  functions  centered  on  each 
molecule.  Their  method  was  based  on  earlier  work  by  Berne  and  Pechukas^^,  and  can  be 
used  to  model  molecules  with  shapes  in  the  form  of  ellipsoids,  both  prolate  and  oblate.  The 
Gaussian  overlap  potential  has  the  advantage  of  being  easily  calculated  and  diflferentiated. 
But  the  CuPc  molecule  cannot  be  well  approximated  as  a  disk-like  ellipsoid;  as  shown  in 
Fig.  1,  it  has  a  nearly  constant  thickness  and  a  definite  four-fold  symmetry  that  we  expect 
plays  an  important  role  in  determining  crystal  structure.  As  an  alternative  to  our  interac¬ 
tion  site  model  one  could  attempt  to  formulate  a  more  complicated  single  site-site  potential 
by  adding  four-fold  symmetric  terms  to  the  Gaussian  overlap  model.  However  such  a  mod¬ 
ification  might  well  compromise  the  model’s  main  strengths  of  being  easily  calculated  and 
differentiated. 

It  is  difficult,  though  not  impossible,  to  construct  a  more  general  single  site-site  potential 
for  a  rigid  molecule  whose  shape  is  arbitrary;  the  Gaussian  overlap  model  is  only  one  example 
of  many  possible  functional  forms.  Another  choice  would  be  to  describe  the  CuPc  molecule 
as,  say,  a  tile-shaped  hard  core  plus  attractive  forces.  However  the  implementation  of  such 
a  potential  in  a  modelling  application  may  prove  difficult  for  other  reasons.  For  instance, 
the  use  of  a  non-spherical  hard  potential  in  a  Monte  Carlo  simulation  requires  a  significant 
amount  of  computer  time  to  test  for  the  presence  of  particle  overlaps  in  trial  moves.  Even 
for  a  soft  potential  written  as  an  explicit  function  of  molecular  separation  and  relative 
orientation,  some  geometry  calculations  are  required  to  determine  the  relevant  angles  and 
distances.  While  we  have  not  fully  explored  the  option  of  constructing  a  single  site-site 
potential,  we  believe  that  the  inherent  simplicity  of  the  interaction-site  reduced  potential 
makes  it  a  better  choice  for  molecules  whose  symmetry  is  neither  spherical  nor  cylindrical. 

Zhang  and  Forrest^^’^"*  followed  yet  another  approach  in  their  work  on  modelling  quasiepi- 
taxial  ordering  at  interfaces  in  organic  molecular  thin  films.  Starting  with  atomistic  po¬ 
tentials  for  PTCDA  (C24O6H8)  and  NTCDA  (CuOeHa)  their  goal  was  to  determine  the 
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degrees  of  freedom.  Atoms  can  form  bonds  that  hold  a  molecule  together  and  if  we  approxi¬ 
mate  the  molecules  as  rigid,  without  internal  degrees  of  freedom,  we  need  not  calculate  these 
intramolecular  interactions.  Assuming  that  there  are  no  chemical  bonds  between  different 
molecules,  the  intermolecular  interaction  consists  of  strong  short  ranged  repulsive  interac¬ 
tions  when  molecular  cores  overlap,  and  at  longer  distances  the  more  slowly  varying  van 
der  Waals  and  electrostatic  interactions.  The  repulsive  and  van  der  Waals  interactions  are 
described  conveniently  together  using  a  Lennard-Jones  potential 


=  4eij 


(2.1) 


where  subscripts  i  and  j  denote  the  species  of  the  atoms  involved  in  the  interaction  and  Vij 
the  distance  between  the  two  atoms.  The  electrostatic  interactions  are  caused  by  the  partial 
charges  of  the  atoms. 

A.  Pohorille*  has  developed  an  atomistic  potential  model  for  CuPc  based  on  the  em¬ 
pirical  energy  functions  and  parameters  used  in  CHARMM.  Henceforth,  we  call  this  the 
CHARMM  model.  The  relative  positions  of  the  atoms  inside  a  molecule  are  determined  by 
X-ray  diffraction  studies^®.  Pohorille  adjusted  partial  charges  within  the  molecule  to  give 
a  potential  which  agrees  well  with  the  /?-CuPc  crystal  structure.  To  calculate  the  inter- 
molecular  interaction  while  neglecting  the  relaxations  of  the  molecules,  we  need  to  sum  over 
all  the  pairs  formed  by  atoms  from  different  molecules;  there  are  57  x  57  =  3249  separate 
interactions  between  two  CuPc  molecules.  As  mentioned  above,  it  is  important  to  have  a 
more  computationally  efficient  model  in  any  crystal  growth  study.  The  approach  we  have 
chosen  is  to  replace  the  atomistic  interaction  by  a  smaller  number  of  “site-site”  interactions. 
For  the  purpose  of  this  exercise  we  will  sissume  the  CHARMM  model  is  exact,  and  try  to 
adjust  the  positions  of  the  sites  and  the  interactions  between  them  so  that  the  intermolec- 
ular  “site-site”  potential  agrees  with  the  CHARMM  model  potential  for  as  many  relevant 
configurations  as  possible. 

For  two  rigid  molecules,  the  configurations  form  a  6-dimensional  continuous  space.  In 
principle  we  could  choose  a  finite  number  of  configurations,  represent  the  positions  and 
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interactions  of  the  sites  by  a  set  of  parameters,  and  use  a  multivariable  minimization  code, 
e.g.  the  simulated  annealing  method^^,  to  find  the  values  of  the  parameters  which  give  the 
best  fit.  However,  the  selection  of  an  “objective  function”,  whose  minimization  corresponds 
to  the  best  fit  of  the  parameters,  is  a  complicated  problem.  We  cannot  fit  the  details  of 
the  short  ranged  repulsion  exactly  in  any  case  and  the  difference  between  the  CHARMM 
model  potential  and  the  “reduced”  model  potential  will  certainly  be  big  for  those  high 
energy  configurations  that  probe  this  highly  repulsive  portion  of  the  interaction.  But  that 
does  not  necessarily  mean  a  bad  fit  in  many  practical  situations  since  those  high  energy 
configurations  contribute  little  to  the  partition  function  or  to  configurations  encountered 
in  crystal  growth.  It  is  more  important  to  fit  those  configurations  which  have  significant 
attraction  correctly  while  reproducing  the  overall  “size”  and  “shape”  of  the  molecule.  We 
do  not  know  any  systematic  way  of  constructing  a  physically  relevant  “objective  function” 
for  so  complex  a  system.  Therefore  our  approach  is  to  first  understand  physically  the  gross 
features  of  interactions  of  the  CHARMM  model.  Then  we  try  to  make  our  “reduced”  model 
refiect  those  features  we  believe  to  be  most  important  by  choosing  parameters  of  the  site-site 
model:  the  positions  and  interactions  between  the  sites. 

The  gross  features  of  the  CHARMM  model  for  CuPc  can  be  separated  into  two  parts: 
1)  the  geometric  “shape”  of  the  molecule  which  is  mainly  determined  by  the  repulsive 
forces,  and  2)  the  attraction  between  molecules.  These  are  closely  related  to  the  prop¬ 
erty  of  Lennard-Jones  potential.  At  short  distances,  the  potential  is  harshly  repulsive  and 
resembles  the  hard  sphere  interaction.  The  electrostatic  and  van  der  Waals  interactions  on 
the  other  hand  are  smoother  at  a  respectable  distance  so  that  the  “size”  and  the  “shape”  of 
the  molecule  is  more  or  less  determined  by  the  repulsive  part  of  the  Lennard-Jones  potential. 
Dij  =  marks  the  distance  where  the  Lennard-Jones  potential  (Eq.  (2.1))  reaches  its 
minimum  and  the  force  changes  from  attractive  to  repulsive.  In  choosing  the  parameters  for 
the  Lennard-Jones  potentials,  the  CHARMM  model  uses  the  combining  rules  which  can  be 
formulated  as 
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^ij  —  i^ii  + 


(2.2) 


€ii  =  (2.3) 

Therefore  in  the  CHARMM  model  we  need  to  know  only  the  interaction  parameters  between 
atoms  of  the  same  species.  Eq.  (2.2)  is  analogous  to  the  additive  rules  for  the  diameters  of 
physical  hard  spheres. 

The  CuPc  molecule  consists  of  57  atoms  which  are  virtually  coplanar  (Fig.  1);  there  are 
seven  chemically  distinct  species  of  atoms.  (A  single  species  of  atoms — e.g.  nitrogen  and 
carbon — can  have  different  properties  depending  on  their  positions  in  the  molecule.)  The 
Da's  for  all  species  are  around  3.7  A  except  for  the  copper  and  hydrogen  atoms  which  are 
about  3.0  A.  Effectively,  we  can  think  of  the  short  ranged  repulsive  part  of  the  CHARMM 
potential  as  made  of  57  partially  overlapping  hard  spheres  of  almost  the  same  size  packed 
together  and  taking  the  “shape”  of  a  square  tile  approximately.  The  tile  has  a  side  length 
about  10  A  and  a  thickness  about  3.5  A.  Of  course  the  molecule  will  be  a  little  bumpy  and 
Fig.  2  illustrate  this  feature  by  representing  each  atom  by  a  sphere  with  diameter 

The  long  ranged  part  of  the  Lennard-Jones  potential  and  electrostatic  interactions  mainly 
contributes  to  the  attraction  between  molecules.  Since  the  CuPc  is  electrically  neutral  and 
has  spatial  inversion  symmetry  about  the  copper  atom,  the  leading  terms  in  the  multipole 
expansion  of  the  local  charge  distribution  are  quadrupoles.  Therefore  asymptotically  the 
interaction  between  two  rigid  and  static  molecules  decreases  as  fast  as  1/r®  at  large  distances. 
The  attraction  is  much  bigger  when  two  molecules  come  face  to  face  than  edge  to  edge. 

To  model  these  features,  we  use  an  assembly  of  “interaction  sites”  and  let  each  of  them 
take  the  role  of  a  cluster  of  atoms  in  the  real  molecule.  Thus  the  atomistic  potential  is 
reduced  to  a  much  smaller  number  of  “site-site”  interactions,  which  we  take  to  be  Lennard- 
Jones  potentials.  In  principle  we  are  free  to  choose  an  arbitrary  function  for  “site-site” 
interaction  but  there  is  little  to  be  gained  by  using  a  more  complicated  function.  Due  to  the 
lack  of  net  electrical  charge  and  dipole  moments  of  an  individual  molecule,  the  Lennard-Jones 
potential  is  physically  reasonable  and  has  an  appropriate  long  ranged  behavior.  Interactions 
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involving  more  than  two-body  terms  would  increase  the  number  of  calculations  tremen¬ 
dously.  Also  by  using  spherically  symmetric  interactions  rather  than  angular  directional 
ones  at  each  site,  we  can  avoid  complications  in  implementation  of  the  reduced  potential  in 
a  modelling  application,  e.g.  molecular  dynamics  study,  by  treating  each  site  as  a  particle. 
The  calculation  of  force  and  torque  will  then  be  straightforward. 

Having  decided  on  the  interaction  sites  approach  and  the  form  of  site  interactions,  the 
next  step  is  to  choose  the  number  and  positions  of  the  sites  and  the  parameters  for  the 
“site-site”  potentials.  As  in  the  CHARMM  model,  we  tried  to  separate  the  fitting  into  two 
parts:  first  to  fit  the  “size”  and  “shape”  of  the  molecule,  then  to  fit  the  attraction  between 
two  molecules.  This  separation  is  important  in  guiding  our  search  of  the  parameters  because 
otherwise  we  would  have  to  go  through  a  multivariable  search  which  is  an  immense  amount 
of  work. 

The  minimum  number  of  the  interaction  sites  is  easily  determined.  The  arrangement  of 
the  interaction  sites  needs  to  satisfy  a  4-fold  rotational  symmetry.  The  Dab's  (we  use  o,  6  •  •  • 
to  denote  the  species  of  “interaction  sites”  with  a  —  A,B,C •  • '  and  Dab's  are  defined  in  the 
same  way  as  those  for  atomistic  interactions)  are  essentially  fixed  by  the  thickness  of  the  real 
molecule,  which  is  about  3.5  A.  Therefore  in  order  to  get  a  side  length  of  about  10  A  without 
leaving  big  “holes”  in  the  middle  we  need  at  least  nine  “interaction  sites” .  The  arrangement 
of  the  sites  (see  Fig.  3(a))  assumes  a  4-fold  rotational  symmetry  and  refiection  symmetry  and 
there  are  three  distinct  species  of  sites.  The  distance  of  each  site  from  the  center  is  so  chosen 
to  give  the  right  size  of  the  molecule.  Thus  in  order  to  provide  an  adequate  description  of 
the  “size”  and  the  “shape”  of  the  molecule,  most  of  the  repulsive  force  parameters  are  fixed 
to  a  certain  range  and  we  have  relatively  little  freedom  to  adjust  them. 

Next  we  fit  the  remaining  parameters,  which  are  essentially  the  eo^’s  of  the  “site-site” 
interaction,  to  the  attractive  part  of  the  CHARMM  model.  Unlike  Eq.  (2.2),  Eq.  (2.3)  is  not 
very  well  substantiated  physically,  and  we  gain  flexibility  by  using  more  general  combining 
rules  for  €ab's  of  the  “interaction  sites” . 

There  are  several  features  of  the  CHARMM  potential  which  deserve  special  attention 
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when  attempting  a  reduction  of  the  potential.  In  Fig.  4(a)  we  show  the  CHARMM  potential 
between  two  parallel  molecules  with  equal  orientation  and  fixed  vertical  displacement  of  3.4 
A,  plotted  as  a  function  of  horizontal  displacement  along  two  coordinates.  It  seems  plausible 
that  such  “sliding”  configurations,  which  maximize  attractive  overlap,  will  be  very  important 
in  the  crystal  growth  process  as  well  as  in  the  statistics  of  equilibrium  crystal  structures. 
Therefore  it  is  very  important  to  fit  these  configurations  as  accurately  as  possible.  For  two 
parallel  molecules  sliding  past  each  other  at  a  fixed  and  not  too  small  distance,  the  interaction 
between  them  is  smooth  in  all  directions.  This  feature  underscores  the  similarity  between 
the  molecule  and  a  square  tile.  The  potential  reaches  its  minimum  when  the  centers  of  the 
two  molecules  are  slightly  displaced.  This  probably  plays  a  role  in  the  frequently  observed 
CuPc  herring-bone  crystal  structure.  Since  for  each  molecule  the  CHARMM  potential  has 
57  interaction  sites,  the  area  where  two  parallel  molecules  have  significant  attractions  is  very 
wide.  In  effect  the  attraction  is  spread  out  over  all  57  atoms  in  the  molecule.  For  our  reduced 
potential  model,  due  to  the  small  number  of  interaction  sites,  it  is  very  hard  to  distribute  the 
attraction  widely  enough,  and  we  could  only  satisfy  this  smoothness  feature  qualitatively. 
Although  it  is  possible  to  fit  any  one  of  these  configurations  much  more  closely,  inevitably 
this  will  reduce  the  quality  of  other  configurations.  Our  aim  is  to  have  a  relatively  robust 
qualitative  model  that  will  work  for  as  many  situations  as  possible.  In  those  studies  that 
don’t  involve  general  configurations,  one  could  adjust  the  parameters  to  fit  those  particular 
configurations  of  interest  better.  To  some  degree  our  model  is  a  compromise  that  includes 
as  much  physics  as  possible  without  using  too  many  sites. 

After  extensive  study  of  this  “9-site-model”,  we  felt  that  there  were  still  significant 
discrepancies  with  the  CHARMM  model,  which  were  especially  severe  for  the  configurations 
when  site-A  of  one  molecule  falls  into  the  gap  between  the  type-B  sites  of  the  other  one. 
We  are  obliged  to  put  four  additional  sites  (see  Fig.  3(b))  to  fill  the  gap  and  since  that 
is  their  only  role,  they  do  not  need  to  interact  with  every  species  of  site;  we  choose  them 
to  interact  with  type-A  sites  only.  That  will  increase  the  number  of  separate  pairs  from 
81  for  the  “9-site-model”  to  89.  The  “13-site-model”  gives  a  potential  which  is  relatively 
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smooth  when  the  molecules  slide  past  each  other  in  all  directions.  Fig.  4(b)  is  the  two 
dimensional  energy  surface  for  the  reduced  model  potential  with  13  interaction  sites.  Prom 
the  spacing  of  the  contours  one  can  see  that  the  potential  minimum  is  much  wider  for  the 
CHARMM  potential  and  we  are  not  able  to  achieve  complete  quantitative  agreement  with 
it.  There  are  simply  not  enough  sites  in  our  model  to  spread  out  the  attractions  enough. 
Fig.  5  shows  the  comparison  of  the  CHARMM  and  reduced  potentials  for  two  representative 
sliding  directions.  Our  model  potential  breaks  down  much  more  severely  when  the  molecules 
are  very  close  to  each  other  and  the  repulsion  is  high  as  illustrated  in  Fig.  6.  However,  since 
these  configurations  are  not  expected  to  appear  often  and  contribute  little  to  the  statistics 
of  the  crystal  structure,  we  think  the  reduced  potential  is  adequate  for  many  studies.  The 
parameters  are  listed  in  Table  I.  The  distance  between  type-A  and  type-S  sites  is  7.0  A, 
and  the  distance  between  type-A  and  type-C  sites  is  4.7  A. 

There  are  three  distinct  stages  in  our  development  of  a  reduced  potential  fit  to  the 
CHARMM  model.  First,  by  choosing  the  right  diameters  and  positions  of  the  sites,  we  can 
get  a  good  description  of  the  “size”  and  "shape”  of  the  molecule.  Then  by  using  appropriate 
“site-site”  interactions,  we  obtain  the  right  long  ranged  behavior.  Finally,  by  fine-tuning 
the  remaining  parameters,  we  achieve  as  good  a  quantitative  agreement  with  the  “exact” 
potential  as  we  can. 


III.  CRYSTAL  STRUCTURES 

Though  we  have  fit  the  gross  features  of  the  reduced  potential  to  the  full  atomistic 
potential,  many  differences  between  the  two  functions  remain.  Here  we  discuss  the  resulting 
differences  in  bulk  crystal  structure,  as  determined  via  energy  minimization. 

To  determine  stable  (and  metastable)  zero  temperature  crystal  structures  for  a  given 
potential,  we  located  minima  in  a  total  energy  function  of  many  variables  which  parametrize 
the  structure  of  the  unit  cell.  We  carried  out  such  a  calculation  for  the  case  of  two  (rigid) 
molecules  per  unit  cell,  for  both  the  CHARMM  potential  and  the  reduced  potential.  The 
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total  number  of  variables  parameterizing  the  crystal  structure  is  18:  nine  corresponding  to 
the  components  of  three  vectors  defining  the  unit  cell;  three  corresponding  to  the  vector 
locating  the  second  particle  in  the  unit  cell;  and  six  corresponding  to  the  orientation  angles 
of  the  two  molecules.  Symmetry  allows  this  number  to  be  reduced.  However  in  hopes  of 
avoiding  barriers  to  finding  local  minima,  we  diose  to  minimize  over  all  18  variables.  This 
choice  did  not  prevent  the  minimization  code^®  from  converging,  though  it  likely  increased 
the  computation  time  unnecessarily. 

For  the  CHARMM  potential,  we  found  two  local  minimum  energy  structures  with  com¬ 
parable  energies.  The  first,  shown  in  Fig.  7(a,b),  is  extremely  close  to  the  known  herringbone 
/3-CuPc  structure  (pictured  in  Ref.^);  this  was  not  unexpected,  as  Pohorille  fit  the  partial 
charges  in  the  CHARMM  model  specifically  to  give  a  low  energy  for  the  /5-CuPc  structure. 
This  herringbone  structure  has  energy=  —120.84  kcal/mole.  We  also  found  a  “flat”  struc¬ 
ture,  shown  in  Fig.  7(c).  Here  the  two  molecules  in  the  unit  cell  are  approximately  degener¬ 
ate,  so  there  is  actually  only  one  molecule/unit  cell.  The  flat  structure  has  energy=— 118.2 
kcal/mole. 

For  the  reduced  potential,  we  also  found  two  structures.  First,  we  found  a  herringbone 
structure,  shown  in  Fig.  8(a,b),  which  resembles — but  is  not  identical  to — the  a-CuPc  struc¬ 
ture  (pictured  in  Ref.^).  This  herringbone  structure  has  energy=-116.1  kcal/mole.  We  also 
found  a  “flat”  structure  in  which  the  two  molecules  in  the  unit  cell  are  approximately  de¬ 
generate,  shown  in  Fig.  8(c).  This  flat  structure  has  energy=— 121.2  kcal/mole.  The  volume 
per  unit  cell  is  somewhat  smaller  for  this  structure  than  for  the  CHARMM  flat  structure 
discussed  above. 

In  working  with  the  reduced  potential,  we  did  not  find  a  local  minimum  close  to  the  j3- 
CuPc  structure.  However  the  search  over  parameter  space  required  in  the  minimization  over 
multiple  variables  is  difficult.  Though  we  started  the  minimization  from  several  different 
initial  conditions — including  the  /3-CuPc  structure  itself — we  cannot  be  sure  that  a  stable 
;9-CuPc-related  crystal  phase  does  not  exist  for  the  reduced  potential.  Likewise,  in  our  work 
with  the  CHARMM  potential,  we  did  not  find  an  a-CuPc  related  crystal,  but  that  may  be 
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because  our  search  was  not  thorough  enough. 

However,  a  closer  look  at  the  exact  form  of  the  two  potentials  yields  further  insight 
into  the  observed  crystal  structures.  In  the  a-CuPc  structure,  molecules  within  a  column 
are  stacked  so  that  the  center  of  each  molecule  is  only  slightly  offset  from  the  one  below; 
while  in  the  ;3-CuPc  structure  the  offset  is  larger.  We  turn  again  to  the  potential  curves 
for  one  molecule  sliding  over  another,  plotted  in  Fig.  5.  From  these  curves  we  see  that  the 
reduced  potential  favors  a  small  offset;  that  is,  energy  is  minimized  when  parallel  molecules 
are  shifted  only  slightly,  as  in  the  a-CuPc  structure.  In  contrast,  the  CHARMM  potential 
favors  a  larger  offset,  as  in  the  /3-CuPc  structure.  While  this  one  aspect  of  local  packing  is 
not  the  only  factor  affecting  the  choice  of  crystal  structure,  it  is  perhaps  the  most  important, 
and  appears  to  explain  why  the  two  potentials  favor  different  structures. 

Roughly  speaking,  crystal  structure  is  determined  by  the  need  to  optimize  attractive 
interactions  while  ensuring  that  the  repulsive  “hard”  cores  do  not  overlap.  In  the  case 
of  CuPc,  these  considerations  favor  stacked  columnar  configurations  where  molecules  are 
packed  as  closely  together  as  the  size  of  the  repulsive  cores  will  allow,  while  maximizing 
attractive  interactions  from  adjacent  molecular  faces.  Since  CuPc  prefers  some  amount 
of  offset  between  molecules  in  a  column,  straight  columns  are  somewhat  less  energetically 
favorable;  columns  can  be  uniformly  tilted,  as  in  the  flat  phases  shown  in  Figs.  6(c)  and  7(c); 
or  the  tilts  can  alternate,  as  in  the  herringbone  phases  shown  in  Figs.  6(a,b)  and  7(a,b).  It 
is  important  to  note  that  these  different  types  of  structures  are  nearly  degenerate  in  energy, 
varying  by  only  a  few  percent,  regardless  of  which  potential  is  studied.  That  is,  as  long  as 
local  packing  is  more  or  less  optimized,  the  rest  of  the  details  of  the  crystal  structure  do 
not  change  the  total  energy  very  much.  This  finding  helps  to  understand  the  experimental 
observation  that  the  crystal  structure  of  a  CuPc  thin  film  grown  from  the  vapor  phase 
depends  sensitively  on  the  structure  of  the  underlying  substrate^^. 

We  conclude  from  this  analysis  that  both  the  CHARMM  and  reduced  potentials  show 
at  least  some  agreement  with  experimentally  observed  crystal  structures.  Both  potentials 
display  herringbone  structures  characteristic  of  Phthalocyanines,  as  well  as  “flat”  phases. 
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Because  x-ray  studies  of  polymorphs  of  CuPc  other  than  at-  and  ;9-CuPc  have  not  been 
performed,  it  is  not  known  if  the  “flat”  structure  corresponds  to  another  known  polymorph^. 

We  doubt  that  any  site-site  potential,  including  atomistic  potentials  such  as  the 
CHARMM  potential,  will  be  able  to  reproduce  consistently  the  small  energy  differences 
that  stabilize  difierent  crystal  polymorphs.  We  conclude  that  the  reduced  potential  ap¬ 
pears  to  be  at  least  roughly  consistent  with  the  type  of  crystal  structure  common  for  CuPc; 
for  this  reason  we  are  optimistic  that  it  may  be  used  in  modelling  applications  to  get  at 
least  approximate  information  about  the  energetics  of  microscopic  processes,  e.g.  surface 
diflFiision. 


IV.  CONCLUSION 

Computer  modeling  of  organic  molecular  thin  films  is  significantly  limited  by  both  com¬ 
puter  processing  time  and  storage  constraints  when  atomistic  potentials  are  used.  In  spite 
of  these  difficulties,  small  systems  have  been  modeled  using  full  atomistic  potentials.  C. 
D.  England  et  al^^  studied  a  model  of  coronene  adsorbed  on  M0S2.  They  note  that  due  to 
limitations  of  both  disk  space  and  processor  capability,  their  computations  were  limited  to 
a  layer  of  seven  coronene  molecules  over  a  substrate  of  size  24  x  24  S-Mo-S  units.  In  related 
work^  the  same  authors  modelled  a  layer  composed  of  9  Pc  molecules  over  a  semiconductor 
substrate.  The  authors  point  out  that  the  small  overlayer  sizes  used  make  it  difficult  to 
draw  any  hard  conclusions  about  the  orientation  angle  and  energies  determined  from  these 
models.  Thus  the  limitations  imposed  by  the  use  of  an  atomistic  potential  are  quite  seri¬ 
ous;  the  use  of  a  reduced  potential  would  dramatically  reduce  both  storage  and  processing 
requirements  and  make  it  possible  to  model  much  larger  system  sizes. 

Using  the  method  described  in  Section  II,  a  computationally  efficient  potential  can  be 
determined  for  any  organic  molecular  material  whose  intramolecular  degrees  of  freedom  can 
be  neglected.  Once  this  task  is  complete,  theory  and  modelling  may  be  used  to  address  a  host 
of  research  questions  concerning  the  growth  and  structure  of  organic  molecular  thin  films. 
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going  beyond  calculations  of  interface  orientations  and  energetics.  As  a  first  step,  crystal 
structures  can  be  determined  through  the  procedure  described  in  Section  III  and  matched 
with  known  polymorphs.  Then,  the  potential  may  be  used  to  calculate  various  properties  of 
the  bulk  crystal,  such  as  elastic  constants,  interlayer  shear  strength,  and  energies/structures 
of  defects.  It  may  also  be  used  to  calculate  static  properties  related  to  the  crystal  surface, 
such  as  anisotropic  surface  energy,  and  energies  of  defects  such  as  steps  or  dislocations.  This 
type  of  calculation  is  especially  important  for  organic  molecular  materials  because  surface 
defects  play  a  critical  role  in  the  response  and  selectivity  of  thin  film  sensors^. 

One  would  also  like  to  implement  computationally  efficient  potentials  to  investigate  the 
microscopic  mechanisms  that  control  crystal  growth.  A  logical  choice  is  to  calculate  the 
energy  barriers  associated  with  diffusion  on  any  given  crystal  face.  This  is  a  very  complicated 
calculation  because  orientational  degrees  of  freedom  as  well  as  center-of-mass  motion  play 
an  important  role  in  diffusion,  such  that  a  molecule’s  effective  diffusion  barrier  depends 
sensitively  on  its  angular  orientation.  We  have  already  begun  such  calculations  for  CuPc  and 
find  that  this  effect  is  very  strong;  this  work  will  be  described  in  a  subsequent  publication^^. 

Use  of  computationally  efficient  potentials  will  also  make  molecular  dynamics  simulation 
feasible.  However,  straightforward  simulations  of  this  type  are  unlikely  to  provide  much 
insight  into  the  crystal  growth  process  because  diffusion  barriers  are  relatively  high,  so  that 
long  time  scales  inaccessible  via  molecular  dynamics  would  be  required  to  observe  significant 
growth.  Instead,  data  from  a  careful  study  of  diffusion  mechanisms  and  barriers  can  be  used 
as  input  parameters  in  a  Monte  Carlo  simulation  in  which  molecules  diffuse  from  one  lattice 
site  to  another  according  to  stochastic  rules,  in  analogy  to  work  done  to  model  kinetic 
formation  of  island  shapes  on  metal  surfaces^.  This  type  of  investigation  is  particularly 
needed  in  light  of  the  recent  experimental  result  of  M.  K.  Debe  and  R.  J.  Poirier^®  showing 
that  uniform  1000  A  thick  films  of  an  organic  red  pigment  will  undergo  a  transition  to  a 
dense  random  array  of  crystalline  whiskers  under  annealing  in  vacuo.  The  authors  assert 
that  in  this  experiment  the  formation  of  whiskers  is  driven  only  by  surface  diffusion,  and 
not  by  kinetic  effects  associated  with  growth,  i.e.  shadowing. 
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With  the  development  of  simplified  potentials  plus  anticipated  gains  in  computer  pro¬ 
cessing  speed,  we  may  soon  be  able  to  simulate  realistic  models  of  these  and  other  such 
complicated  processes  in  organic  molecular  thin  films.  It  appears  likely  that  work  in  this 
area  will  lead  to  new  understanding  of  the  fundamental  properties  of  these  materials. 
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FIGURES 


FIG.  1.  Chemical  constituents  of  CuPc. 

FIG.  2.  Picture  of  the  CuPc  molecule  representing  each  atom  by  a  sphere. 

FIG.  3.  Positions  of  the  “interaction  sites”  for  (a)  9-site  model  and  (b)  13-site  model. 

FIG.  4.  The  two  dimensional  dimer  energy  for  (a)  CHARMM  model  and  (b)  reduced  potential 
model  for  two  parallel  molecules  with  the  same  orientation. 

FIG.  5.  Comparison  of  CHARMM  model  (solid  line)  and  reduced  potential  model  (dashed  line) 
for  two  parallel  molecules  with  interplanar  distance  of  3.4  A.  (a)  shows  the  comparison  when  the 
molecules  slide  along  the  diagonal  direction  of  Fig.  1  or  3;  (b)  x-axis. 

FIG.  6.  Comparison  of  CHARMM  model  (solid  line)  and  reduced  potential  model  (dashed  line) 
for  two  parallel  molecules  with  interplanar  distance  of  2.8  A. 

FIG.  7.  Crystal  structxires  obtained  by  finding  local  minima  of  total  intermolecular  energy 
function  of  a  unit  cell.  For  CHARMM  potential,  two  difierent  structmres  are  foimd:  (a)  and  (b) 
show  the  first  structure,  which  is  close  to  the  herringbone  /3-CuPc  structure,  from  two  viewing 
directions;  another  “fiat”  structure  which  is  shown  in  (c)  is  also  foimd. 

FIG.  8.  (a)  and  (b)  show  the  herringbone  structure  found  by  energy  minimization  using  the 
reduced  potential.  It  is  close  to — ^but  not  identical  to — the  a-CuPc  structme.  As  in  the  case  of 
CHARMM  model,  a  “flat”  structure  (c)  is  also  found  for  the  reduced  potential,  where  the  two 
molecules  in  the  unit  cell  axe  approximately  degenerate. 
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TABLES 


TABLE  1.  Parameters  for  the  site-site  interactions  for  the  reduced  potential  model. 


intereiction  type 

D{k) 

e{kcal/mol) 

A- A 

4.0 

-0.8 

A-B 

3.6 

-1.8 

A-C 

3.6 

-1.9 

B-B 

3.8 

-0.5 

B-C 

3.6 

-1.0 

C-C 

3.42 

-3.0 

A-D 

3.65 

-2.0 
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